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

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

CC BY
158
17
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ИСПАРИТЕЛЬ / ХЛАДАГЕНТ / ВНУТРИТРУБНОЕ КИПЕНИЕ / ЭНЕРГОЭФФЕКТИВНОСТЬ / ОПТИМИЗАЦИЯ / ЧИСЛЕННЫЙ МЕТОД / ОГРАНИЧЕНИЕ ПАРАМЕТРОВ / EVAPORATOR / REFRIGERANT / IN-TUBE BOILING / ENERGY EFFICIENCY / OPTIMIZATION / NUMERICAL METHOD / PARAMETER LIMITATION

Аннотация научной статьи по физике, автор научной работы — Кошелев Сергей Валерьевич, Ейдеюс Альгирдас Иозапович

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

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

Похожие темы научных работ по физике , автор научной работы — Кошелев Сергей Валерьевич, Ейдеюс Альгирдас Иозапович

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

ACCOUNT OF LIMITS DURING OPTIMIZATION OF REFRIGERANTS IN-TUBE BOILING PARAMETERS IN SHIPBOARD EVAPORATORS

High heat transfer rate and low pressure drop for a given boiling conditions of the refrigerant is achieved when the difference between the mean temperature of the inner wall of the pipes tТ and the saturation temperature t02 at the outlet from the boiling zone becomes minimal. It takes into account both the irreversible losses caused by the presence of a finite temperature difference, and the effect of the evaporator’s hydraulic resistance on the performance of the compressor. Existing optimization techniques for this criterion only apply to traditional refrigerants and use obsolete formulas to determine the heat transfer coefficient. The numerical method allows to take into account more factors influencing the process of boiling both traditional and alternative refrigerants. Using the computer program in a rather wide range of initial data determining the boiling conditions of the four refrigerants, optimal values of the mass velocity (wρ)о are chosen at known heat flux densities q. After the modification of the program, the optimum length of boiling zone values lко are chosen for the given values of the heat load Qз and other combinations of the initial data. By regression analysis of the chosen results, power dependences are obtained for an approximate determination of the optimal parameters. The values of the power-law coefficients for calculating the basic variables sought are presented in two tables of text. Based on the analysis of the recommendations published in the domestic literature, formulated following the results of an experimental and analytical study of efficiently operating evaporators, it was proposed to introduce, when optimizing, restrictions on the temperature difference tТ and t02, and also on the vapor velocity leaving the boiling zone. To remain within the permissible limits, it is recommended to change the values q and Qз, as well as the internal diameter of the pipes dT and the vapor conten x1 at the refrigerant inlet to the evaporator.

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

ВЕСТНИК«)

ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА ^^

МОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА

DOI: 10.21821/2309-5180-2017-9-4-795-805

ACCOUNT OF LIMITS DURING OPTIMIZATION OF REFRIGERANTS IN-TUBE BOILING PARAMETERS IN SHIPBOARD EVAPORATORS

S. V. Koshelev, A. I. Eideyus

Kaliningrad State Technical University, Kaliningrad, Russian Federation

High heat transfer rate and low pressure drop for a given boiling conditions of the refrigerant is achieved when the difference between the mean temperature of the inner wall of the pipes tT and the saturation temperature t02 at the outlet from the boiling zone becomes minimal. It takes into account both the irreversible losses caused by the presence of a finite temperature difference, and the effect of the evaporator's hydraulic resistance on the performance of the compressor. Existing optimization techniques for this criterion only apply to traditional refrigerants and use obsolete formulas to determine the heat transfer coefficient. The numerical method allows to take into account more factors influencing the process of boiling both traditional and alternative refrigerants. Using the computer program in a rather wide range of initial data determining the boiling conditions of the four refrigerants, optimal values of the mass velocity (wp)o are chosen at known heat flux densities q. After the modification of the program, the optimum length of boiling zone values lo are chosen for the given values of the heat load Q3 and other combinations of the initial data. By regression analysis of the chosen results, power dependences are obtained for an approximate determination of the optimal parameters. The values of the power-law coefficients for calculating the basic variables sought are presented in two tables of text. Based on the analysis of the recommendations published in the domestic literature, formulatedfollowing the results of an experimental and analytical study of efficiently operating evaporators, it was proposed to introduce, when optimizing, restrictions on the temperature difference tT and t02, and also on the vapor velocity leaving the boiling zone. To remain within the permissible limits, it is recommended to change the values q and Q, as well as the internal diameter of the pipes dT and the vapor conten Xj at the refrigerant inlet to the evaporator.

Keywords: evaporator, refrigerant, in-tube boiling, energy efficiency, optimization, numerical method, parameter limitation.

For citation:

Koshelev, Sergey V., and Algirdas I. Eideyus. "Account of limits during optimization of refrigerants in-tube

boiling parameters in shipboard evaporators." Vestnik Gosudarstvennogo universiteta morskogo i rechnogo

flota imeni admirala S. O. Makarova 9.4 (2017): 795-805. DOI: 10.21821/2309-5180-2017-9-4-795-805.

УДК 621.565

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

С. В. Кошелев, А. И. Ейдеюс

ФГБОУ ВО «Калининградский государственный технический университет»,

Калининград, Россия

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

9

4

ГТёТ

ЛВЕСТНИК

............ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Х^ОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА

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

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

Для цитирования:

Кошелев С. В. Учет ограничений при оптимизации параметров внутритрубного кипения хладагентов в судовых испарителях / С. В. Кошелев, А. И. Ейдеюс // Вестник Государственного университета мор -ского и речного флота имени адмирала С. О. Макарова. — 2017. — Т. 9. — № 4. — С. 795-805. DOI: 10.21821/2309-5180-2017-9-4-795-805.

Введение

Энергоэффективность холодильных систем постоянно находится в поле зрения специалистов. На нее влияют многие факторы, среди которых затронем лишь вид хладагента и условия его кипения в испарителе. Как отмечается в статье [1], после принятия Монреальского протокола 1987 г. турбулентности с хладагентами в индустрии холода возникают беспрерывно. Применяется целый ряд озонобезопасных хладагентов, хотя ни один из них не удовлетворяет всем предъявляемым требованиям. Испарители с внутритрубным кипением обеспечивают малую емкость хладагента и относительно высокую энергоэффективность, но их проектирование требует повышенной точности расчета термогидродинамических характеристик [2]. Ввиду сложности процессов кипения жидкостей в трубах и каналах продолжается их исследование и разработка методов прогнозирования режимов течения и истинных параметров двухфазных потоков [3], [4]. Совершенствование методов расчета способствует более точному определению оптимальных параметров кипения хладагентов в испарителях.

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

Разность между средними температурами внутренней стенки труб tт и кипящего хладагента зависит от плотности теплового потока q и среднего КТО аа. Ее можно выразить как ^ - t0 = q/аа. Падение давления в трубах сопровождается понижением температуры насыщения хладагента от ^ на входе до t02 на выходе из зоны кипения. Понижение температуры насыщения Ats = t0l - Оптимальный режим работы достигается при минимуме обеих разностей. Чтобы объединить их, в качестве критерия оптимальности используется разность между температурой стенки tт и температурой хладагента на выходе из зоны кипения

К - = V Ч + у • Ч. (!)

Множитель у учитывает характер понижения температуры кипения по ходу движения хладагента. В статье [5] принято у = 0,5. Авторы статей [6], [7] этот множитель определяют по формуле

ВЕСТНИК,

ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

МОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА,

у = (3-Ах)/(6 - 3Ах), (2)

где Ах = х2 — х1 — приращение массового паросодержания хладагента в зоне кипения.

В статье [8] отмечается, что для большинства практических случаев подходит множитель у = 0,6. При малых значениях Ats = t0l — t02 вполне допустимо принимать у = 0,5. Актуальность нового подхода к оптимизации обусловлена недостатками существующих методик и неизбежным переходом на применение альтернативных хладагентов.

Исходными данными для оптимизации массовой скорости хладагента в трубах ^р являются расчётная температура кипения t0, массовое паросодержание на входе х1 и на выходе х2 зоны кипения, внутренний диаметр труб dт и плотность теплового потока д. Часто внутритрубное кипение хладагента происходит в плоских змеевиках, сформированных из горизонтальных труб, последовательно соединённых калачами. В таких случаях к исходным данным относятся также число труб в змеевике пт, радиус осевой линии калачей Rк и угол наклона плоскости змеевика 0. После подбора оптимальной скорости (^р)о нетрудно определить длину зоны кипения /к и тепловую ее нагрузку Qз. На практике нередко удобнее считать известной тепловую нагрузку Qз при заданных или выбранных значениях х1, х2, dт, а также пт, Rк и 0 для змеевика. В таком случае оптимизация параметров внутритрубного кипения хладагента сводится к подбору длины зоны кипения / , при которой достигается минимум разности ^ - t02 Значения массовой скорости ^р и плотности теплового потока д становятся искомыми переменными.

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

Методы

Опубликованные в статьях [5] - [7] методики оптимизации направлены на поиск оптимального значения массовой скорости (^р)о в трубах испарителей по минимуму разности ^ - t02, при определении которой несколько различается множитель у в формуле (1). Более существенные различия обусловлены тем, что авторы использовали неодинаковые уравнения для определения КТО и падения давления при кипении хладагентов R12 и R22 в горизонтальных трубах. За счёт ряда упрощений каждому автору удалось получить свое уравнение для расчета оптимального значения массовой скорости (^р)о. Недостаток указанных методик обусловлен тем, что их нельзя распространить на новые хладагенты, переход на использование которых вытекает из Международных соглашений по охране окружающей среды.

Методика [8] предполагает оптимизацию длины зоны кипения /к и падения давления хладагента АР. Для определения КТО в ней используется уравнение, которое учитывает изменение массового паросодержания хладагента Ах. Предлагается использовать общий коэффициент сопротивления трению, учитывающий трение в трубах, калачах, а также ускорение потока. Порядок его определения не обсуждается. С учетом характера изменения КТО аа и падения давления АР в зависимости от массовой скорости ^р получено соотношение у •Ats = 0,27 • q / аа. Из него при у = 0,6 получается, что понижение температуры насыщения должно составлять Ats < 0,45 • q / аа. С этим соотношением не согласуется приведенный в статье пример и формулировка, согласно которой понижение температуры Аts должно быть в четыре раза меньше разности температур tт - t0 = д /аа. Окончательное выражение для оптимальной длины / с учетом введенных ранее обозначений может быть записано в виде

-1(1/2,9)

2 О

7

СО

к

4

I = 0,56

а °'8 г3'4

Х(" - V"7°

¿Т-29 Ах °'83

■ а0,62( у¥\°,344

(3)

где ц — динамическая вязкость жидкости; 1 — коэффициент теплопроводности жидкости; г = С - V— удельная теплота парообразования хладагента; V' и V" — удельный объем жидкости

и пара; Т0 — абсолютная температура кипения; Г — общий коэффициент трения, обычно равный 0,02 - 0,03.

Массовая скорость ур не входит в последнее выражение, хотя используется для расчета числа Рейнольдса — Re. Нетрудно заметить, что при прочих равных условиях длина / пропорциональна д-0,62. в данной методике, как и в методиках [5] - [7], по умолчанию предполагается турбулентный режим течения кипящей жидкости, хотя он вполне может быть переходным и даже ламинарным. При определенных условиях происходит расслоение жидкой и паровой фаз хладагента. Ввести условные переходы в одно итоговое выражение практически невозможно. Поэтому предлагается перейти к численному методу подбора оптимальной скорости (ур)о или длины / .

Исследованию гидродинамики и теплообмена при кипении разных хладагентов в трубах посвящено много работ. В последние годы предпочтение отдается обобщенным методикам расчета локальных КТО и градиентов давления. На их основе составлена компьютерная программа [9]. Она распространяется на кипение десяти хладагентов в плоских змеевиках с числом горизонтальных труб от 2 до 20. Основные особенности программы кратко освещены в статьях [10], [11]. Исходными данными для неё служат вид хладагента, температура кипения паросодержания х1 и х2, диаметр труб число труб пт, радиус калачей Rк, угол их наклона ©, плотность теплового потока д и массовая скорость ур. Основными результатами расчета являются средний КТО аа и полное падение давления в змеевике ДРп, хотя попутно определяются длина /к, тепловая нагрузка Q , скорость пара на выходе из зоны кипения уп, понижение температуры насыщения Д^ и ряд других параметров.

Чтобы подбирать оптимальную скорость (ур)о при заданных условиях кипения хладагента, в окно пользователя программы выведены два значения разности ¿т - ¿02. Одно из них получается по формуле (1) при у = 0,5, а другое — вычислением у по формуле (2). Программа [9] определяет полное падение давления ДРп как сумму отдельных составляющих на коротких участках труб змеевика. Тем самым учитываются слабо выраженная криволинейность понижения давления Р0 и температуры ¿0 кипения. Ввиду ограничений на понижение температуры Д^ представляется приемлемым множитель у = 0,5. До перехода к подбору оптимальной длины /ко приведем основные соотношения между показателями кипения хладагента в трубах:

/К = мр • <г (х2 - X )/(4д) = ^ / (); (4)

4Qз 4д • /к ™р—2-з-=-:-• (5)

г (х2 - х1) dT г (х2 - х1)"

Qз= qпdт • /к = • г (х2 - х1 )/4 . (6)

Отсюда следует, что при известных условиях кипения хладагента, используя выбранный критерий оптимальности, можно по заданной плотности теплового потока д подбирать скорость (ур)о и находить соответствующие ей значения /к, Qз, массового расхода хладагента Оа, объемного расхода образующегося пара Vn и скорости пара на выходе из зоны кипения уп. Если задана тепловая нагрузка змеевика Qз, то можно подбирать оптимальную длину /ко и находить соответствующие ей значения д, ур,Оа, V, Уп и других переменных.

Для определения длины /ко по заданной нагрузке Qз составлена модификация программы [9]. К исходным данным в ней вместо д и ур отнесены Qз и длина горизонтальных труб змеевика /т. Искомыми переменными наряду с прежними величинами становятся д и ур, а длина зоны кипения определяется как / = /тпт. Эту модификацию для краткости обозначим как /». По обеим разновидностям программы необходимо последовательно изменять варьируемую величину, т. е. ур при заданной д, или /т при заданной Qз, до достижения минимальной разности ¿т - ^02. Пока этот процесс не автоматизирован. Ввиду плавного приближения к минимуму он оказывается трудоемким, несмотря на почти мгновенное получение результата после ввода очередной переменной. Для ускорения процесса использована программа «Поиск решения». Опыт показывает, что результат по ней нередко зависит от заданного значения варьируемой переменной.

ВЕСТНИК«)

ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА ^^

МОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА

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

Результаты

С использованием программы [9] выполнен подбор оптимальной скорости (wp)o для хладагентов R404A, R410A, R290 и R22 при разных сочетаниях исходных данных: q, dт, t0, xl Они охватывают диапазон q = 500 ... 10000 Вт/м2, dт = 0,005 ... 0,02 м, t0 = -45 ... -5 оС, x1 = 0,05 ... 0,4. Во всех вариантах принято x2 = 0,99, пт = 10, радиус калачей задавался как R = 2 dr., рассматривалось расположение змеевика в вертикальной плоскости при нижней подаче жидкости. Количество неповторяющихся сочетаний исходных данных для каждого хладагента превышало 260. При каждом сочетании наряду с (wp)o получены значения aа, APn, /к, Ats, i, - t02, wn и ряда других переменных. Обобщение многочисленных результатов подбора проводилось с использованием подпрограммы классической линейной модели множественной регрессии. Наиболее подходящей для аппроксимации оказалась степенная зависимость вида

(wp)0 = 10«° q d?(t° + 90)«3 x. (7)

Путем логарифмирования и замены переменных она приводится к линейному уравнению регрессии. Встроенная в пакет Excel функция ЛИНЕЙН после ввода массива исходных данных выдает статистические оценки коэффициентов а. и некоторые показатели точности аппроксимации. Среди них выделим коэффициент детерминированности г2, который должен стремиться к единице, и стандартную ошибку для оценки искомой величины sey, которая должна стремиться к нулю. Степенная зависимость (7) подходит и для аппроксимации других переменных, соответствующих оптимальной скорости (wp)o, но с другими значениями коэффициентов а..

В табл. 1 приводятся коэффициенты уравнения (7) лишь для расчёта (wp)o, /к, аа и tT - t02. При их использовании надо следить, чтобы значения исходных данных не выходили за указанные пределы. Остальные искомые переменные, включая Ats и Q, можно найти по соотношениям (4) и (6) с использованием данных о термодинамических свойствах хладагента на линии насыщения. Скорость выходящего из зоны кипения пара wn = wp • v". Нетрудно определить и полное падение давления APn по найденному понижению температуры Ats относительно расчётной t0. Заметим, что остальные методики оптимизации не выдают в явном виде значений аа и t - t02, хотя полученные по ним данные о q и wp позволяют найти /к и Q.

Таблица 1

Основные итоги регрессионного анализа результатов подбора скорости (wp)

Переменные Статистические оценки

а0 а1 а2 а3 а4 r2 se y

R404

wp -0,56033 0,33165 0,109574 0,927951 0,026291 0,92768 0,042902

/ к 4,253288 -0,66726 1,142993 0,735723 -0,17933 0,976179 0,046473

a а -0,02093 0,542639 0,084079 0,837992 0,087445 0,931246 0,058906

К - t02 1,345793 0,354338 -0,03111 -1,22078 -0,05682 0,935844 0,042449

R410

wp -0,93179 0,357174 0,044701 1,006704 0,021959 0,965928 0,030247

/к 4,121621 -0,63812 1,074081 0,734097 -0,20083 0,987266 0,032384

a а -0,02341 0,512887 -0,06207 0,793528 0,076397 0,947805 0,047534

t, - t02 1,221179 0,378193 0,037432 -1,17086 -0,04703 0,957423 0,036391

2 1

7

ГИ9

ЛВЕСТНИК

............ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Х^ОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА

Таблица 1 (Окончание)

Я290

ур -1,18068 0,384067 0,16381 1,039124 -0,02255 0,915844 0,054698

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

/ к 3,824598 -0,61323 1,19172 0,914251 -0,24426 0,957119 0,060096

а а -0,25559 0,578975 0,095453 0,920975 0,053118 0,924292 0,068004

К - ^>2 1,271888 0,34658 -0,00309 -1,16049 -0,06471 0,925559 0,044615

Я22

ур -1,21615 0,377104 0,105878 1,155761 0,013268 0,927965 0,051315

/ к 3,612476 -0,62343 1,102152 0,970541 -0,17944 0,964763 0,054663

а а -1,26714 0,640161 0,048011 1,271581 0,0791 0,940879 0,070378

К - ^2 1,980877 0,30207 -0,02286 -1,42142 -0,05546 0,957369 0,034394

Подбор оптимальной длины /ко проводился с использованием модифицированной программы /т». После ввода в неё тепловой нагрузки Qз и остальных параметров, определяющих условия кипения выбранного хладагента, подбиралась длина прямых труб /то до достижения минимальной разности ^ - t02. Длина зоны кипения /ко = /то • пТ. Такой подбор выполнен для хладагентов R404A, R410A, R507 и R22. По каждому из них рассмотрено по 350 неповторяющихся сочетаний исходных данных в диапазоне ^ = -45 ... -5 оС; = 0,004 ... 0,02 м; х1 = 0,01 ... 0,4; Qз = 200 ... 8000 Вт. По-прежнему, х2 = 0,99; Як = 2 • dт; пт= 10; угол наклона плоскости змеевика © = 0о; подача жидкости в нижнюю трубу змеевика.

Расчёты показывают, что значения разностей Дts, ^ - существенно зависят от диаметра dт и нагрузки Qз. Чтобы получить приемлемые результаты с уменьшением диаметра dт и понижением температуры следует уменьшать нагрузку Qз. Основными результатами расчета наряду с / являются значения а , ДР , Дt, t - t, д, ур, у . Для компактного представления взаимосвязи

ко а7 п' ^ т а' ' ' п ~ г- ^

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

/ = 10Ьо Qb dbт2 х^(Г0 + 90)ь.

(8)

Такая зависимость оказалась подходящей для аппроксимации не только длины /ко, но и соответствующих ей значений других искомых переменных. Понятно, что для каждой переменной при кипении любого хладагента получаются свои значения коэффициентов Ь . В табл. 2 приводятся коэффициенты Ь0... Ь4 для расчёта лишь значений /ко, аа, ур и ¿т - t02 при кипении четырех хладагентов. Эти данные наряду с нагрузкой Qз позволяют найти д и ур по соотношению (6). Значения коэффициента аа представляют самостоятельный интерес. Они позволяют найти среднюю разность ^ - t0 при оптимальной длине /ко.

Таблица 2

Основные итоги регрессионного анализа результатов подбора длины I

Переменные

Статистические оценки

Ъ.

Ъ.

Ъ

г2

не

Я404

7,739485

-1,17562

3 452007

-0 34029

2,184998

0,959215

0,062498

ур

-5,2546

0,994037

-1,97044

0,163248

0,219576

0,99024

0,023873

а

-6 73887

1482514

-2,91207

0,291571

-0 19918

0,941921

0 083014

К - ^2

—1,27196

0,752669

-1,62158

0,075293

-2,1158

0,960475

0,037398

Ъ

Ъ

г

/

Таблица 2 (Окончание)

Я410

/ ко 7,942196 -1,16517 3,453571 -0,34464 2,161622 0,966707 0,056904

^р -5,36921 0,992145 -1,96791 0,163097 0,21904 0,990083 0,024046

а а -6,88646 1,490796 -2,97351 0,296652 -0,21709 0,955482 0,072959

К - ^ -1,26023 0,698798 -1,50263 0,061957 -1,99541 0,959044 0,036167

Я507

/ ко 7,700551 -1,13681 3,381965 -0,36417 2,05781 0,952935 0,065896

^р -5,25571 0,994663 -1,97098 0,16336 0,223584 0,990308 0,023876

а а -6,68962 1,459958 -2,88172 0,320998 -0,14355 0,938906 0,085003

-1,36113 0,755588 -1,62824 0,077074 -2,08434 0,959385 0,037727

Я22

/ ко 7,841163 -1,23134 3,566147 -0,35578 2,32884 0,959402 0,06461

^р -5,22331 0,991386 -1,97186 0,1679 0,162782 0,990763 0,023004

а а -7,40713 1,59948 -3,17953 0,316035 -0,34812 0,945055 0,086568

К - -0,78201 0,718299 -1,52709 0,072298 -2,17635 0,952318 0,041585

Достаточно высокие значения коэффициента детерминированности г2 при малой ошибке se , в табл. 1 и 2 указывают на хорошую степень корреляции между исходными данными и искомыми переменными, если исходные данные находятся в пределах своего диапазона. Сочетания исходных данных, при которых разности Ats и ^ - заметно превышали 5 °С, не входят в число вариантов, рассмотренных при подборе как скорости (^р) о, так и длины /ко.

Примеры определения оптимальной скорости хладонов R22 и R12 приводятся в статьях [5], [6]. Представляет интерес сравнение приведённых в них данных с результатами численного подбора ^р. В статье [5] недостаточно полно указаны условия кипения хладонов. Более полную информацию ее автор привел в источнике [12]. Для сопоставления по приведенным в источниках [12] и [6] рисункам найдены значения оптимальной скорости хладона R22 при одинаковых плотностях теплового потока q = 2000; 5000; 10000 Вт/м2 и полном испарении жидкости, т. е. х2 = 1. Наряду с данными х1 и они приводятся в табл. 3. Рядом с ними указаны значения (^р)о, подобранные численным методом при тех же условиях кипения. Дополнительно приведены соответствующие им значения длины /к, коэффициента аа и разности tт -

По данным табл. 3 нетрудно заметить, что при t0 = 0 оС приведенные в источниках [12], [6] оптимальные скорости ^р несколько превышают их значения, подобранные численным методом. Предположительно причиной такого превышения может служить приближённый подход к определению потерь давления двухфазного потока в змеевиковых испарителях. В методиках [5], [6] есть лишь упоминание о калачах, а падение давления в них отдельно не рассматривается. Правда, используется фиксированный коэффициент сопротивления трению, который должен учитывать все составляющие падения давления потока.

В случае кипения хладона R22 при t0 = -30 оС оптимальные скорости ^р по данным [12] г

оказались заметно ниже значений, подобранных с использованием программы [9]. Заметим, что по обеим методикам с понижением температуры t0 уменьшается оптимальные скорости ^р. Более резкое уменьшение скорости ^р может быть связано с использованием в книге [8] формул для определения среднего коэффициента теплоотдачи аа, которые в явном виде не учитывают температуру t0 и изменение паросодержания х2 - х1. Расхождение результатов расчета коэффициента аа по двум формулам в источнике [12] объясняется трудно учитываемыми факторами, что свидетельствует об ориентировочном характере полученных данных об оптимальной скорости.

ю

4

ЛВЕСТНИК

............ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Х^ОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА

Таблица 3

Оптимальная скорость R22 по трем методикам

Данные источников [12], [6] Численный подбор

Источник t ос 0' ^ Х1 ^ м д, Вт/м2 ур, кг/(м2с) ур, кг/(м2с) /, м к' аа, Вт/(м2 • К) t -1 ос \ 02' ^

2000 140 132,2 48,9 1680 1,79

[12] 0 0,15 0,017 5000 202 172,5 25,5 2194 2,80

10000 280 205,8 18,3 3014 4,07

2000 57 80,6 38,8 921,2 3,37

[12] -30 0 0,017 5000 90 115,9 22,3 1679 4,38

10000 122 140,5 13,5 2402 5,46

2000 125 113,8 24,8 1475 1,81

[6] 0 0,15 0,01 5000 190 170,8 14,9 2260 2,79

10000 240 240,2 10,5 3337 3,77

2000 140 120,8 52,6 1396 1,88

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

[6] 0 0,15 0,02 5000 217 180,8 31,5 2178 2,87

10000 280 253,4 22,1 2954 4,16

Оптимальная длина труб змеевика /ко определена в статье [12] при следующих условиях кипения хладона R22: = 0 оС; х1 = 0,25; х2 = 1; dт = 0,011 м; д = 5000 Вт/м2. Приняты коэффициент трения F = 0,02 и множитель у = 0,6. По формуле (3) получена длина / = 14,4 м. Попутно найдена скорость пара на выходе из зоны кипения = 8 м/с. Используя программу [9], при тех же исходных данных находим оптимальную скорость (ур) = 171 кг/(м2с). Ей соответствуют искомые переменные /к = 14,44 м; = 8,05 м/с; аа = 2301,4 Вт/(м2К); АРп = 19950 Па; А ^ = 1,25 оС; - = 2,7976 оС; Qз = 2494 Вт. Подставляя значения 2з = 2494 Вт и /т = 1,444 м при пт= 10 в модифицированную программу /т», получаем прежние значения искомых переменных. Это подтверждает правильность модифицированной программы и свидетельствует о хорошем совпадении результатов подбора (ур)о по программе [9] с данными статьи [8].

Если по модифицированной программе /т» проводить оптимизацию длины /к, варьируя длиной /т, то при прежних условиях кипения и Qз = 2494 Вт после поиска решения имеем длину / = 31,54 м, которой соответствуют значения аа = 2141,8 Вт/(м2К), АРп = 32530 Па, д = 2289,4 Вт/ м2, ур = 171 кг/(м2с), Аts = 2,45 оС, tт - t02 = 2,2939 оС. Как видим, при той же нагрузке Qз длина /ко увеличивается в 2,18 раз, а тепловой поток д уменьшается во столько же раз. Несмотря на увеличение АРп в 1,63 раза и А^ в 1,96 раз, показатель оптимальности tт - уменьшается в 1,22 раза. Полученные данные оказались несколько неожиданными. Поэтому были выполнены сравнительные расчеты при других сочетаниях исходных данных. Выяснилось, что во всех случаях подбор (ур)о по заданной д дает более низкие значения /к и АРп, чем при подборе /ко по нагрузке Qз, найденной в результате подбор (ур)о. Расхождения между этими и другими сравниваемыми величинами уменьшаются по мере снижения плотности теплового потока д. Характерно, что значения ур и уп совпадают во всех сравниваемых вариантах.

Отмеченные расхождения обусловлены ограничениями на изменение разных параметров в двух используемых программах. Программа [9] предполагает неизменность плотности теплового потока д, а минимум разности ^ - находится путем подбора массовой скорости (ур)о. По -путно заметим, что на практике постоянство плотности теплового потока д обеспечивается лишь при электрообогреве труб. В остальных случаях она определяется условиями теплопередачи. В модифицированной программе /т» заданной считается тепловая нагрузка Qз, а величина д при подборе длины /ко может изменяться. Минимальная разность (1) при совпадающих значениях Qз, ур и оказывается несколько меньше при использовании программы /т» за счет увеличения длины /к. Анализ сравнительных расчётов показывает, что при подборе длины / уменьшается составляющая д/а уравнения (1), но увеличивается величина Аt по сравне-

нию с результатами подбора (^р)о при аналогичных условиях кипения и одинаковых значениях Q . Выходит, что при подборе / растет коэффициент обратимости обратного цикла [13], но снижается производительность компрессора [14].

Обсуждение

Как следует из ранее изложенного, минимум разности ^ - t02 может быть достигнут при разных сочетаниях параметров на стороне хладагента в зависимости от вида задаваемой величины, т. е. q или Qз. Без учета интенсивности теплоотдачи на стороне охлаждаемой среды значения этих величин можно предсказать лишь приближенно. Обе они зависят от вида охлаждаемой среды, скорости ее движения, а также типа и степени наружного оребрения труб. По-видимому, при проектировании змеевиковых испарителей можно по аналогии с существующими задавать q и находить (^р)о.

Подбор типоразмера змеевиков из числа выпускаемых промышленностью применительно к конкретным условиям кипения хладагента лучше проводить в зависимости от нагрузки Qз. В общем случае, целесообразно по ориентировочным значениям q найти (^р)о и Qз, а затем по Qз определить / . На основе сопоставления результатов можно принять окончательное решение с учетом других ограничений. Уже отмечалось, что нежелательны сочетания исходных данных, при которых значения ^ - t02 и/или Д^ превышают 5 оС. Как правило, в таких случаях скорость насыщенного пара wn превышает 15 м/с. Таким значением, в частности, ограничивают скорость пара в линии всасывания компрессора. В зарубежной практике эту скорость ограничивают так, чтобы понижение температуры насыщения на участке от испарителя до компрессора не превышало 2^ = 1,1 оС. С позиций влияния на производительность компрессора подобные ограничения следует распространить и на зону кипения хладагента в трубах. Исходя из возможных описок в статье [8], представляется целесообразным упомянутое ранее ограничение записать как Дts < 0,45(^ - ^2). Для рассмотренного в ней примера по программе [9] получено Дts = 1,25 оС, а 0,45(^ - ^2) = 1,26 оС, что вполне правдоподобно.

По экспериментальным и расчетным данным в отечественной литературе приводятся рекомендуемые значения средней логарифмической разности ©т между температурами охлаждаемой среды и кипящего хладагента. В большинстве случаев она находится в пределах ©т = 6 ... 10 оС при охлаждении как жидкости, так и воздуха. С понижением температуры t0 разность ©т несколько уменьшается. Рассмотренные испарители чаще всего работают при массовой скорости хладагента wр = 40 ... 140 кг/(м2с). Для судовых холодильных установок нормальной считается температура кипения t0, которая на 5 - 7 оС ниже средней температуры жидкого хладоносителя и на 7 - 10 оС ниже температуры охлаждаемого воздуха. В автономных кондиционерах, зная температуру охлаждаемого воздуха tъl, принимают температуру кипения как t0 = tвl + (13 ... 17 оС), но добиваются средней температуры наружной поверхности ребристого воздухоохладителя tн = t0 + (5 ... 7 оС). Поскольку температурный напор ©т имеет три составляющие, предположим, что одну треть его может занимать разность tт - t0 = (6...10 оС)/ 3. Тогда q/аа = 2 ... 3,3 оС. Не должна превышать этого значения и величина Дts, обусловленная падением давления Д Р Если считать, что проведенную в статье [8] устную формулировку следует записать как Дt = © /4, то получается Дt < 1,5 ... 2,5 оС.

Вытекающие из разных рекомендаций ограничения не совпадают. Полагая допустимыми значения q /аа = 2 ... 3,3 оС и Дts 1,5 ... 2,5 оС, следует ограничить разность ^ - в пределах 2,75 ... 4,55 оС, а скорость пара wn — 15 м/с. С понижением температуры t0 из-за снижения плотности пара р" при прочих равных условиях увеличиваются значения ДРп, Дt и wn. По -этому оптимальная скорость ^р)о уменьшается. Тем не менее, как следует из табл. 3, при t0 = -30 оС разность ^ - становится выше 4,55 оС, когда q — 5000 Вт/м2. Чтобы учесть обсуждаемые ограничения, приходится уменьшать задаваемы значения q или Qз, хотя их увеличение желательно для интенсификации теплоотдачи на стороне хладагента и уменьшения числа змеевиков в испарителе.

2 О

7

СО

г

4

ЛВЕСТНИК

............ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Х^ОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА

Заключение

Оптимизация параметров внутритрубного кипения хладагентов по минимуму разности ^ - t02 может привести к неверным результатам, если неправильно выбрать исходные данные. Предлагается добиваться значений ^ - t02 < 4,55 оС, ее составляющих д/аа < 3,3 оС, Аts < 2,5 оС и скорости насыщенного пара уп < 15 м/с путем ограничения плотности теплового потока д или тепловой нагрузки на змеевик Qз, а также подбора диаметра труб dт с учетом температуры кипения t0 и начального паросодержания хг Приведенные численные значения нуждаются в дальнейшем уточнении; требуется подготовка рекомендаций по определению оптимального числа змеевиков, направленных на минимизацию массогабаритных показателей испарителя, что актуально для судового оборудования.

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

Составленная программа [9], в отличие от существующих методик [5] - [8], наряду с оптимальным значением массовой скорости (ур)о или длины зоны кипения / выдает значения остальных показателей процесса кипения до десяти хладагентов. Путем регрессионного анализа результатов подбора (ур)о и / для четырех хладагентов получены степенные зависимости, позволяющие найти как (ур)о или / , так и соответствующие им основные показатели. Значения коэффициентов степенной зависимости приводятся в табл. 1 и 2 для использования в инженерной практике. Планируется дополнить программу модулем автоматического поиска минимальный разности ^ - t02 и выдачи предупреждения о выходе искомых переменных за установленные пределы.

СПИСОК ЛИТЕРАТУРЫ

1. Цветков О. Б. Энерго- и экологически эффективные рабочие вещества в технологиях генерации холода и теплоты / О. Б. Цветков, Ю. А. Лаптев // Холодильная техника. — 2016. — № 3. — С. 18-24.

2. Малышев А. А. Перспективные типы испарителей холодильных машин / А. А. Малышев, В. О. Мам-ченко, В. М. Мизин [и др.] // Вестник Международной академии холода. — 2013. — № 2. — С. 13-18.

3. Малышев А. А. Новые методы прогнозирования режимов течения кипящих хладагентов в макро- и миниканалах / А. А. Малышев, К. В. Киссер, А. С. Филатов // Вестник Международной академии холода. — 2016. — № 2. — С. 67-70.

4. Малышев А. А. Истинные параметры кипящих хладагентов в трубах и каналах / А. А. Малышев, К. В. Киссер, А. В. Зайцев // Вестник Международной академии холода. — 2017. — № 2. — С. 53-56.

5. Гоголин А. А. Об оптимальной скорости фреона в трубках испарителей / А. А. Гоголин // Холодильная техника. — 1965. — № 1. — С. 29-33.

6. Захаров Ю. В. Определение оптимальной массовой скорости хладагента в горизонтальных трубках испарителей / Ю. В. Захаров, Н. И. Радченко // Холодильная техника. — 1980. — № 3. — С. 25-29.

7. Granryd E. Om val av seriekopplad rörlängd och tryckfall vid förongare med fullstandig förongning / E. Granryd // Kylteknisk tidskrift. — 1966. — № 4. — Pp. 65-68.

8. Granryd E. Optimum Circuit Tube Length and pressure drop on the Refrigerant side of evaporator / E. Granryd // International Refrigeration and Air Conditioning Conference. — Paper 143. — 1992. — Pp. 73-82.

9. Свидетельство о государственной регистрации программы для ЭВМ №2015663262. Общая программа расчета коэффициента теплоотдачи и падения давления при кипении десяти хладагентов в плоских змеевиках с разным числом труб / А. И. Ейдеюс, С. В. Кошелев; правообладатель ФГБОУ КГТУ. — Дата регистрации 31.08.2015 г.

10. Ейдеюс А. И. Сравнение интенсивности теплоотдачи и падения давления при кипении хладагентов R404A и R22 в горизонтальных трубах / А. И. Ейдеюс, М. Ю. Никишин, C. В. Кошелев // Вестник Международной академии холода. — 2015. — № 1. — С. 69-74.

11. Ейдеюс А. И. Теплоотдача и падение давления при кипении хладагентов в змеевиках / А. И. Ейдеюс, C. В. Кошелев, М. Ю. Никишин // Вестник Международной академии холода. — 2016. — № 2. — С. 42-47.

12. Гоголин А. А. Интенсификация теплообмена в испарителях холодильных машин/ А. А. Гоголин, Г. Н. Данилова, В. М. Азарсков, Н. М. Медникова. — М.: Легкая и пищевая промышленность, 1982. — 224 с.

ВЕСТН1

ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

МОРСКОГО И РЕЧНОГО ФЛОТА ИМЕНИ АДМИРАЛА С. О. МАКАРОВА

13. Кошелев С. В. Влияние условий внутритрубного кипения хладагентов на внешнюю необратимость в испарителях / С. В. Кошелев, А. И. Ейдеюс, М. Ю. Никишин // Известия КГТУ. — 2016. — N° 42. — С. 117-125.

14. Ейдеюс А. И. Падение давления в змеевиковых испарителях и производительность компрессора / А. И. Ейдеюс, С. В. Кошелев, М. Ю. Никишин // Компрессорная техника и пневматика. — 2015. — № 5. — С. 14.

1. Tsvetkov, O. B., and Yu. A. Laptev. "Energoi ekologicheski effektivnye rabochie veshchestva v tekhnologi-yakh generatsii kholoda i teploty." Kholodil'naya tekhnika 3 (2016): 18-24.

2. Malyshev, Aleksandr A., Valery O. Mamchenko, V. M. Mizin, A. V. Potanina, S. I. Proshin, and T. I. Devya-tov. "Perspective types of evaporators of refrigerators." Vestnik Mezhdunarodnoi Akademii Kholoda 2 (2013): 13-18.

3. Malyshev, A. A., K. V. Kisser, and A. S. Filatov. "New methods of forecasting flow regimes for boiling refrigerant in macro and mini channels." Vestnik Mezhdunarodnoi Akademii Kholoda 2 (2016): 67-70.

4. Malyshev, A. A., K. V. Kisser, and A. V. Zaitsev. "true parameters of boiling refrigerants in tubes and channels." Vestnik Mezhdunarodnoi Akademii Kholoda 2 (2017): 53-56.

5. Gogolin, A. A. "Ob optimal'noi skorosti freona v trubkakh isparitelei." Kholodil'naya tekhnika 1 (1965): 29-33.

6. Zakharov, Yu. V., and N. I. Radchenko. "Opredelenie optimal'noi massovoi skorosti khladagenta v gorizontal'nykh trubkakh isparitelei." Kholodil'naya tekhnika 3 (1980): 25-29.

7. Granryd, E. "Om val av seriekopplad rorlangd och tryckfall vid forongare med fullstandig forongning." Kylteknisk tidskrift 4 (1966): 65-68.

8. Granryd, E. "Optimum Circuit Tube Length and pressure drop on the Refrigerant side of evaporator." International Refrigeration and Air Conditioning Conference. Paper 143. 1992. 73-82.

9. Eideyus, A. I., and S. V. Koshelev. Certificate for computer programm №«2015663262. Obshchaya programma rascheta koeffitsienta teplootdachi i padeniya davleniya pri kipenii desyati khladagentov v ploskikh zmeevikakh s raznym chislom trub. Russian Federation assignee. Publ. 31 Aug. 2015.

10. Eideyus, Algirdas I., Mikhail Yu. Nikishin, and Sergey V. Koshelev. "The comparison of the intensity of heat transfer and pressure drop during boiling of the refrigerants R404A and R22 in horizontal pipes." Vestnik Mezhdunarodnoi Akademii Kholoda 1 (2015): 69-74.

11. Eideyus, A. I., S. V. Koshelev, and M. Yu. Nikishin. "Heat transfer and pressure drop at refrigerant boiling in the plate coils." Vestnik Mezhdunarodnoi Akademii Kholoda 2 (2016): 42-47.

12. Gogolin, A. A., G. N. Danilova, V. M. Azarskov, and N. M. Mednikova. Intensifikatsiya teploobmena v isparitelyakh kholodil'nykh mashin. M.: Legkaya i pishch. prom-st', 1982.

13. Koshelev, C. V., A. I. Eideyus, and M. Yu. Nikishin. "Vliyanie uslovii vnutritrubnogo kipeniya khladagentov na vneshnyuyu neobratimost' v isparitelyakh." Izvestiya KGTU 42 (2016): 117-125.

14. Eideyus, A. I., S. V. Koshelev, and M. Y. Nikishin. "Pressure Drop in Coil Evaporators and Compressor Capacity." Kompressornaya tekhnika ipnevmatika 5 (2015): 14.

REFERENCES

ИНФОРМАЦИЯ ОБ АВТОРАХ

INFORMATION ABOUT THE AUTHORS

Кошелев Сергей Валерьевич — аспирант

Научный руководитель:

Ейдеюс Альгирдас Иозапович

ФГБОУ ВО «Калининградский государственный

технический университет»

236022, Российская Федерация, Калининград,

Советский проспект, 1

e-mail: entermoria@rambler.ru

Ейдеюс Альгирдас Иозапович —

кандидат технических наук, профессор

ФГБОУ ВО «Калининградский государственный

технический университет»

236022, Российская Федерация, Калининград,

Советский проспект, 1

e-mail: xktk@bga.gazinter.net

Koshelev, Sergey V. — Postgraduate

Supervisor:

Eideyus, Algirdas I.

Kaliningrad State Technical University

1 Sovetskij Prospect, Kaliningrad, 236022,

Russian Federation

e-mail: entermoria@rambler.ru

Eideyus, Algirdas I. —

PhD, professor

Kaliningrad State Technical University 1 Sovetskij Prospect, Kaliningrad, 236022, Russian Federation e-mail: xktk@bga.gazinter.net

Статья поступила в редакцию 13 июля 2017 г.

Received: July 13, 2017.

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