Полный текст

Введение В настоящее время проблема безопасной утилизации отработанного ядерного топлива (ОЯТ) все больше выходит на передний план. Ежегодные объемы, извлекаемые из реакторов только в Российской Федерации, составляют около 700 т. В США объем накопленного ОЯТ составляет порядка 40 тысяч т. При сохранении современных тенденций проблема безопасного хранения и накопления ОЯТ неминуемо приведет к техногенной катастрофе. Технологии переработки ОЯТ существуют, но они дороги и практически нивелируют преимущества ядерной энергетики. Единственным рентабельным способом переработки ОЯТ является производство топлива для реакторов на быстрых нейтронах. Разработка технологии замкнутого ядерного цикла ведется достаточно давно и практически близка к завершению. Значимой частью данной технологии является пирохимический передел. Данная операция заключается в высокотемпературной обработке нитридного СНУП-топлива в составе опытно-демонстрационного энергетического комплекса (ОДЭК). Пирохимическая часть - одна из наиболее технологически сложных и инновационных составляющих проекта замкнутого ядерного цикла. Основная задача пирохимического передела - сокращение времени выдержки ОЯТ до начала переработки, чтобы как можно быстрее запустить ценные компоненты в ядерный топливный цикл, а также сократить общее количество пребывающих в нем делящихся материалов [1]. Реализация пирохимического процесса накладывает жесткие ограничения по изолированности среды, где находится продукт переработки. Для этого используются инертные камеры (ИК). ИК, как правило, устроена как замкнутая двустенная оболочка, во внутреннем и межстеночном объеме которой находится аргон. При этом в межстеночном пространстве поддерживается давление ниже атмосферного. Тем самым при нарушении герметичности как внутренней, так и наружной оболочек исключается проникновение рабочей среды в атмосферу [2]. Основное требование к ИК, помимо герметичности, - обеспечение съема тепла, выделяемого обрабатываемым продуктом и обеспечение температур оболочек и технологического оборудования внутри ИК в допускаемых пределах. Так как рабочий процесс ИК не подразумевает нахождение человека во внутреннем объеме в течение всего жизненного цикла оборудования, все операции, как технологические, так и ремонтные, осуществляются дистанционно управляемым манипулятором. Следовательно, все компоненты манипулятора, такие как уплотнения, система смазки и т.д., должны выдерживать температурные режимы в течение всего срока эксплуатации. Основная задача исследования - разработка и валидация численной теплофизической модели, которая позволит прогнозировать тепловое состояние проектируемых ИК и принимать конструктивные решения для обеспечения требуемых температур ограждающих конструкций камеры и технологического оборудования в ней при максимизации количества обрабатываемого продукта. Материалы и методы Для разработки теплофизической модели использовалась CAD-модель ИК, представленная на рис. 1. Ограждающая конструкция ИК выполнена в виде двух оболочек из нержавеющей стали и силового каркаса между ними из конструкционной стали обыкновенного качества. Система охлаждения представлена в виде пяти блоков охлаждения, размещенных на внутренних стенах. Данные блоки являются центробежными вентиляторами, которые прокачивают внутреннюю среду камеры через теплообменник (рис. 2). Рис. 1. CAD-модель ИК (разрез продольный) а б Рис. 2. Блок охлаждения (а) и роботизированный манипулятор (б) Привод крыльчаток блоков охлаждения осуществляется извне, через магнитную муфту. Хладагент прокачивается через теплообменник отдельной системой, размещенной снаружи ИК. Самым массивным и объемным оборудованием, находящимся во внутреннем объеме ИК, является роботизированный манипулятор (см. рис. 2), перемещающийся по направляющим, размещенным в верхней части камеры. Манипулятор обеспечивает выполнение всех требуемых технологических операций: перемещение контейнеров с продуктом, открытие реторт, извлечение реторт и складирование на стеллажах и т.д. Также манипулятор обеспечивает все необходимые ремонтные процедуры: замена подшипников крыльчаток блоков охлаждения, замена камер видеонаблюдения и элементов внутреннего освещения и т.д. В качестве источника тепла, имитирующего тепловыделение продукта, использовался массогабаритный тепловой макет (рис. 3). Для нагрева применялись 16 электрических тэнов общей мощностью 10 кВт с возможностью плавной регулировки. Для контроля температуры во внутреннем объеме и межстеночном пространстве ИК были размещены 15 датчиков (рис. 4). Рис. 3. 3D-модель макета массогабаритного теплового а б Рис. 4. Расположение основных датчиков на продольной стене камеры (а), на торцовых стенах (б) Часть дополнительных датчиков размещалась на вертикальных струнах во внутреннем объеме камеры (рис. 5). Также под потолком на уровне опорных кронштейнов, направляющих перемещения манипулятора, расположены еще два дополнительных датчика (рис. 6). Рис. 5. Расположение дополнительных датчиков Рис. 6. Расположение дополнительных потолочных датчиков Время выполнения экспериментального исследования составляло 25 ч, температура в цеху, где размещалась инертная камера, поддерживалась на уровне 23,2 °С, максимальная мощность нагревателей макета - 10 кВт, среда во внутреннем объеме и ограждающей конструкции камеры - воздух. Мощность, подаваемая на нагревательные элементы макета, регулировалась в диапазоне от 0 до 100 % от максимальной мощности. Измерения мощности проводились для четырех групп нагревателей, по четыре нагревателя в группе. График изменения мощности макета приведен на рис. 7. Как видно из рис. 7, фактически наблюдалась некоторая неравномерность потребляемой мощности по группам нагревателей. Данные по показаниям основных и дополнительных датчиков до момента отключения теплового макета (38 000 с) представлены на рис. 8. Рис. 7. График изменения мощности теплового макета а б Рис. 8. Показания датчиков основных (а) и дополнительных (б) Как видно из графиков, ИК не была доведена до равновесного теплового состояния. Нагрев ИК продолжался вплоть до выключения массогабаритного теплового макета. Данное обстоятельство несколько снизит точность валидации численной модели и определит необходимость моделирования нестационарного процесса, который требует значительных расчетных и временных ресурсов. Для снижения вычислительной сложности расчетов детализация полости межстеночного пространства ИК, в которой проходит инертная среда, и технологическое оборудование внутри камеры упрощены до простых геометрических форм. Смотровые окна и технологические проемы имеют двойные стенки, между которыми поддерживается избыточное давление инертной среды, в связи с этим смотровые окна и технологические проемы выделены в отдельные замкнутые объемы в ячейках. Исходная геометрия оболочки, силового каркаса и фрагмент расчетной геометрии представлен на рис. 9. Рис. 9. Исходная геометрия оболочки, силового каркаса и фрагмент расчетной геометрии Общий вид упрощенной геометрии ограждающей конструкции камеры представлен на рис. 10. Рис. 10. Упрощенная геометрия ограждающих конструкций камеры инертной В состав технологического оборудования, учитываемого в расчетной модели, вошли: манипулятор силовой и путь рельсовый, макет массогабаритный тепловой, блоки охлаждения (рис. 11). Рис. 11. Расчетная геометрия внутреннеего объема КИ с манипулятором, блоками охлаждения и тепловым макетом В результате упрощения исходной геометрии получена расчетная геометрия оболочки, внутреннего объема и основного технологического оборудования, участвующего и/или влияющего на тепловое состояние. Толщины стенок элементов силового каркаса и оболочек, их теплоемкость и теплопроводность будут учитываться в расчетной модели через назначение специальных интерфейсов с функцией Shell Conduction. Такой подход позволит адекватно учесть теплопередачу вдоль стенок тонкостенных элементов конструкции. Методика моделирования КИ базируется на следующих допущениях: - теплопроводность и теплоемкость материалов, применяемых для ограждающих конструкций камеры инертной, не меняются во всем рабочем диапазоне температур; - условия окружающей среды неизменны в течение всего времени выполнения испытаний; - потери электрической мощности в кабелях от измеряющего и регулирующего устройства до массогабаритного теплового макета отсутствуют. Этапы выполнения расчетов следующие: а) выполнение расчетов, обосновывающих пригодность упрощенной геометрии для применения при расчетах камеры инертной; б) выполнение расчетов фрагмента ограждающей конструкции камеры для определения возможности выполнения нестационарного расчета напрямую либо необходимости введения допущений, снижающих ресурсоемкость моделирования; в) выполнение расчета внутреннего объема с учетом ограждающих конструкций в виде граничных условий, обоснованных по результатам выполнения п. «б»). Основная идея упрощенной геометрии - исключение из расчетной модели оболочек и силового каркаса как твердых тел и замена их поверхностями нулевой толщины. При этом математически толщина, теплопроводность и теплоемкость в расчете учитываются в виде одномерных моделей. Применение функции Shell Conduction также позволяет учитывать тепловые потоки не только по нормали к стенке, но и вдоль ее плоскости. Моделирование теплового режима камеры инертной проводилось в программе Ansys Fluent. Во внутреннем пространстве камер рассматривалось однофазное движение сжимаемого идеального газа. Учитывался теплообмен за счет конвекции, теплопроводности и излучения. Процессы в газообразной среде описываются основными уравнениями классической гидродинамики [3]. Для учета турбулентного характера движения газа в КИ использовалась двухпараметрическая модель турбулентности k-ε [4], являющаяся на сегодня инженерным стандартом. В модели турбулентности использовались стандартные значения констант и настроек [5]. Учет лучистого теплообмена осуществлялся с помощью модели Surface-to-Surface (S2S). Модель предполагает оптическую прозрачность среды, то есть учитывает только излучение от поверхности к поверхности и игнорирует поглощение, испускание и рассеяние лучистой энергии в среде между ними. Все поверхности полагаются диффузными и серыми. Расчет коэффициентов видимости производился с помощью алгоритма трассировки лучей с количеством лучей, равным 50, кластеризация не использовалась [6]. Математическая модель замыкается заданием граничных условий на поверхностях расчетной области. В расчетах задавались следующие динамические и тепловые граничные условия: · условие заданного массового расхода газа (вход газа); · условие заданного статического избыточного давления (выход газа); · условие прилипания (все поверхности кроме входа и выхода газа); · условие заданной температуры (поверхности оборудования, вход газа); · условие смешанного (конвекция и излучение) теплообмена с окружающим пространством с заданной температурой (внешняя поверхность камер). Условие прилипания задавалось на всех границах газового объема камер (стенках камеры и поверхности оборудования). Для источника тепловой энергии в составе макета задавалась мощность тепловыделения на поверхности нагревателей. Для внешних поверхностей камеры внешняя температура и коэффициент теплоотдачи принимается в соответствии с результатами вспомогательного расчета, исходя из температуры окружающей среды = 23,2 °С. Так как степень черноты оказывает значительное влияние на параметры лучистого теплообмена [7], согласно общепринятому инженерному опыту приняты следующие значения данного параметра: на стенках камеры = 0,4; на внешних поверхностях оборудования = 0,8; на поверхностях внутренних полостей оборудования = 1. Расчеты проводились с использованием сопряженного решателя по давлению. Полная расчетная сетка ИК строилась при помощи сеточного генератора Fluent meshing. Построение велось по ярусам с контролем правильности построения с точки зрения стенок и зазоров. Внешний вид и поперечное сечение построенной расчетной сетки представлены на рис. 12. Рис. 12. Внешний вид расчетной сетки ограждающей конструкции Количество элементов расчетной сетки составило чуть меньше 20 млн элементов. Внутренний объем ИК (см. рис. 11) моделировался отдельно. Внешний вид расчетной сетки представлен на рис. 13. Рис. 13. Расчетная сетка внутреннего объема камеры Размер расчетной сетки внутреннего объема камеры составил 3,25 млн элементов. С учетом типа сетки - polyhedral, можно констатировать тот факт, что расчетные мощности, требуемые для выполнения нестационарного расчета длительностью 38 000 с - до момента отключения теплового макета (см. рис. 8, б), - будут значительны, и время, требуемое на выполнение данного расчета, будет исчисляться месяцами. Исходя из получаемых размеров расчетных сеток, прямая постановка задачи затруднительна с точки зрения требуемых ресурсов и требуемого времени для выполнения нестационарных расчетов. Единственный способ моделирования теплового состояния ИК в данной ситуации - отказ от моделирования ограждающей конструкции и приведение его к граничным условиям на внутренней оболочке камеры. Для определения отказа от рассмотрения ограждающей конструкции определим число Био для оболочек ограждающей конструкции ИК: где - коэффициент конвективной теплоотдачи, Вт/(м2·К); - толщина оболочки; - теплопроводность для материала оболочек. Для определения коэффициентов конвективной теплоотдачи и тепловых потоков были выполнены предварительные расчеты углового фрагмента ограждающей конструкции ИК (см. рис. 9), визуализация результатов представлена на рис. 14. а б в Рис. 14. Поля распределения температур при различных температурах внутренней стенки: а - 80 °С; б - 51,6 °С; в - 57,5 °С Средние расчетные коэффициенты конвективной теплоотдачи для крыши, наружных стен и пола представлены в табл. 1. Таблица 1 Средние расчетные коэффициенты конвективной теплоотдачи, Вт/(м2К) Параметр Значение Крыша 11,22 Стенка 6,5 Дно 4,2 Результаты расчета числа Био приведены в табл. 2. Таблица 2 Число Био для элементов наружной оболочки Параметр Значение Крыша 0,0016 Стенка 0,0028 Дно 0,001 В соответствии со значениями числа Био оболочки ограждающей конструкции камеры можно рассматривать как тонкие тела и не рассматривать градиент температур на их поверхностях. Удельный расчетный тепловой поток через внешнюю оболочку представлен в табл. 3. Таблица 3 Удельный тепловой поток через элементы наружной оболочки, Вт/м2 Параметр Значение Крыша 76,3 Стенка 50 Дно 5 Далее определим термическое сопротивление для боковых стен, потолка и пола по формуле где - средняя температура внутренней оболочки, °С; - средняя температура внешней оболочки, °С; - удельный расчетный тепловой поток, Вт/м2. Средняя температура зон наружной оболочки представлена в табл. 4. Таблица 4 Средняя температура зон наружной оболочки, °С Параметр Значение Крыша 37 Стенка 34,3 Дно 25,1 Результаты расчета термического сопротивления приведены в табл. 5. Таблица 5 Термическое сопротивление элементов ограждающей конструкции, К/Вт Параметр Значение Крыша 0,269 Стенка 0,464 Дно 6,48 Отсюда получим приведенные коэффициенты теплопроводности для определения граничных условий в расчетной модели по формуле где - толщина ограждающей конструкции по средней линии оболочек. Результаты расчета приведенных коэффициентов теплопроводности представлены в табл. 6. Таблица 6 Приведенные коэффициенты теплопроводности, Вт/(м·К) Параметр Значение Крыша 0,789 Стенка 0,457 Дно 0,033 Далее получим приведенные коэффициенты конвективной теплоотдачи для определения граничных условий в расчетной модели по формуле Результаты расчета приведенных коэффициентов конвективной теплоотдачи представлены в табл. 7. Таблица 7 Приведенные коэффициенты конвективной теплоотдачи, Вт/(м2К) Параметр Значение Крыша 4,488 Стенка 3,497 Дно 0,98 Таким образом, дальнейшие расчеты выполнялись только с использованием расчетной сетки внутреннего объема (см. рис. 13). Процессы теплообмена с окружающей средой через ограждающие конструкции камеры будут учтены посредством наложения граничных условий в виде приведенных коэффициентов теплоотдачи, приведенных значений теплопроводности и соответствующей толщины стенок ограждающей конструкции камеры через назначение интерфейсов. Теплопередача вдоль элементов ограждающей конструкции будет осуществляться путем иcпользования функции Shell conduction. Как показали предварительные расчеты, при соответствии расчетных температур реторты массогабаритного макета экспериментальным мощность нагревательных элементов макета составляет 1,83 кВт. При этом в эксперименте данные температуры достигались при полной мощности макета, равной 4 кВт. Данное обстоятельство связано с тем, что в реальной конструкции макета на корпусе имеются отверстия (см. рис. 3), через которые происходит конвективный съем тепла из внутреннего объема макета. Площадь нижних отверстий составила 4396 мм2, а верхних - 7435 мм2. Так как данный процесс смоделировать достаточно проблематично, то учтем данную часть съема тепла с макета путем итерационного увеличения теплопроводности теплоизоляционного материала макета: ТКП-1350 (теплопроводность от 0,23 до 0,3 Вт/(м·К)). Увеличение теплопроводности позволит перераспределить температуры по объему макета и увеличить конвективную теплоотдачу с внешних поверхностей макета. Данный подход не нарушает теплового баланса всей камеры и не повлияет на точность выполнения расчетов. Для этого было выполнено проектирование модели макета, состоящего из трех частей: реторта, зона нагревательных элементов и корпус. Модель и расчетная сетка показаны на рис. 15. Рис. 15. Модель теплового макета Таким образом, была определена методика и численные модели для расчета теплового состояния ИК и их валидации на основе полученных экспериментальных данных. Результаты Стационарный расчет проводился по параметрам условно стационарного режима испытаний макета, который длился с 30 000 с до момента отключения теплового макета - 38 000 с (см. рис. 7). Согласно результатам предварительных расчетов теплового макета, итерационно было определено, что теплопроводность теплоизоляции макета, равная 0,78 Вт/(м·К), обеспечивает значение температуры в районе нагревателя 600 °С, что при мощности макета 4 кВт соответствует экспериментальным данным. Стационарный расчет проводился до момента установления температуры всех датчиков. Результаты расчета и сравнение их с экспериментальными данными представлено в табл. 8. Таблица 8 Расчетные и экспериментальные значения датчиков температуры Обозначение Значение Отклонение, % расчет протокол расчет протокол ts1-low ТС-1 33,1 25,82 22,1 ts2 ТС-2 35,7 36,01 -1,0 ts3 ТС-3 35,0 36,2 -3,5 number1-hight №1 40,7 40,4 0,6 number2-middle №2 36,2 34,5 4,7 number3-low №3 36,0 30,2 16,1 number4-roof №4 42,5 41,1 3,3 number5-roof №5 41,4 40,1 3,1 number6-hight №6 40,8 40,5 0,8 number7-hi/mid №7 40,3 40,3 -0,1 Значительное отклонение датчика ТС-1 объясняется тем, что в эксперименте он был закрыт от излучения макета инструментальным стеллажом, который был перенесен из другого места внутреннего объема камеры (рис. 18). Отклонение дополнительного датчика № 3 объясняется тем фактом, что для учета конвективной теплоотдачи внутри макета была изменена теплопроводность теплоизоляции макета, что привело к повышению температуры стенок макета. А так как датчик № 3 находится в непосредственной близости от стенки макета (см. рис. 5) - это привело к его дополнительному нагреву от излучения. Рис. 16. Изменение положения инструментального стеллажа Рис. 17. Поле распределения температур Рис. 18. Поле распределения скоростей конвективных потоков В целом расчетные значения датчиков несколько выше экспериментальных, так как режим работы камеры в период времени от 30 000 до 38 000 с не является стационарным, и разогревание конструкции продолжалось. Распределение температур скоростей конвективных потоков в продольном сечении показано на рис. 17 и 18. Нагрев внутренней оболочки потолка показан на рис. 19. Рис. 19. Поле распределения температур на внутренней оболочке потолка камеры Полученные данные по температурам внутренней оболочки камеры свидетельствуют, что температура на наружной оболочке не достигнет 40 °С. Далее был выполнен нестационарный расчет для валидации расчетной модели с точки зрения динамики нагрева ИК. В расчете учитывалась динамика изменения мощности теплового макета от времени (см. рис. 7). Результаты расчета в виде графиков нагрева датчиков представлены на рис. 20. а б Рис. 20. Расчетные и экспериментальные данные по основным датчикам Как видно из графиков, расчетные значения датчиков хорошо согласуются с данными эксперимента. Расхождение практически по всем датчикам не превышает 5 %, за исключением начального этапа моделирования (до 6000 с). Значительное расхождение значений пар датчиков ТС-1, ts1-low и Т-3, number3-low объясняется теми же причинами, что и в стационарном расчете. Для предварительного определения эффективности системы охлаждения ИК и максимальной производительности по количеству тепловых макетов был выполнен расчет теплового состояния камеры с тремя макетами и работающей системой охлаждения. Соответствующая геометрия внутреннего объема приведен на (рис. 21). Рис. 21. Расчетная область с тремя массогабаритными тепловыми макетами В конструкции камеры предусмотрены пять охлаждающих блоков. Так как данных по зависимости температуры на выходе в блок охлаждения от температуры на входе нет, была принята температура воздуха на выходе из блока охлаждения - 10 °С, перекачиваемый объем аргона - 30 м3/ч. Результаты расчета в виде значений температур датчиков представлены в табл. 9. Таблица 9 Расчетные и экспериментальные значения датчиков температуры Обозначение датчика Значение, °С Расчет Протокол ts1 ТС-1 24,2 ts2 ТС-2 19,1 ts3 ТС-3 20,5 number1-hight №1 32,9 number2-middle №2 28,3 number3-low №3 26,1 number4-roof №4 20,4 number5-roof №5 21,9 number6-hight №6 25,3 number7-hi/mid №7 27,6 Поле распределения температур на внутренней оболочке потолка камеры представлено на рис. 22, а. а б в Рис. 22. Поле распределения температур: а - на внутренней оболочке потолка камеры; б - боковых стенок и пола камеры; в - боковых стенок со стороны блоков охлаждения Максимальная расчетная температура потолка камеры не превышает 30 °С. Нагрев боковых стенок и пола показан на рис. 22, б. Продольная стенка, противоположная от стены с блоками охлаждения, является самой теплой, максимальная температура - 31,5 °С. Соответственно, стенки со стороны блоков охлаждения является самой холодной (рис. 22, в), их температура не превышает 28 °С. Скорости потоков охлажденной среды не превышают 3,7 м/с (рис. 23). Распределение температур по поверхности манипулятора представлено на рис. 24. Рис. 23. Поле скоростей потоков охлажденной среды Рис. 24. Поле распределения температур на поверхности манипулятора Наиболее нагретая часть манипулятора находится практически над открытой ретортой 3-го макета. Смещение макетов в противоположную сторону позволит снизить темповую нагрузку на колонну манипулятора. Также, если положение манипулятора постоянно, можно направить поток охлажденной среды от блока охлаждения, находящегося на торцевой стенке, на колонну манипулятора. Общее распределение температур в продольном сечении представлено на рис. 25. Рис. 25. Распределение температур в продольном сечении Результаты расчета показывают, что система охлаждения с приведенными параметрами обеспечивает требуемые параметры теплового состояния камеры инертной. Обсуждение Разработанная численная теплофизическая модель ИК для пирохимического передела ОЯТ обладает достаточной адекватностью для применения при разработке и проектировании данного оборудования. Отличительной характеристикой данной модели являются граничные условия, моделирующие ограждающую конструкцию без ее прямого одновременного моделирования. Данное обстоятельство позволяет моделировать как стационарное, так и нестационарное тепловое состояние ИК с использованием меньших вычислительных мощностей. При прямом моделировании полной ограждающей конструкции ИК получается очень большой размер конечно-элементной сетки, так как ограждающая конструкция имеет малые конструктивные элементы (зазоры между силовым каркасом и оболочками составляют 4 мм), которые необходимо учитывать при моделировании. Поэтому граничные условия разработанной модели обосновывались посредством моделирования части фрагмента ограждающей конструкции и определения ее теплотехнических характеристик для дальнейшего применения в качестве граничных условий при моделировании теплового состояния ИК в целом. Для более точного моделирования тепловых состояний пирохимического оборудования необходима полная информация о комплектующих: расходно-напорные характеристики, зависимости интенсивности теплообмена в зависимости от расходов и температур как хладагента, так и охлаждаемой среды и т.д. При выполнении экспериментальных исследований для более точной валидации численных теплофизических моделей необходимо достигать стационарных состояний. Так как значительная часть тепловой энергии передается от тепловых источников излучением, необходима полная информация о степени черноты поверхностей элементов конструкции пирохимического оборудования. Заключение Таким образом, результатом данной работы является валидированная методика математического моделирования теплового состояния двухоболочечных ИК и оборудования пирохимического передела в части применимости допущений и упрощений математической модели. Выполненные экспериментальные исследования динамики нагрева исследуемой ИК с использованием массогабаритного теплового макета показали хорошую согласованность экспериментальных и расчетных данных. Разработанная методика может быть применима для дальнейшего проектирования оборудования пирохимического передела ОЯТ, что позволит обеспечить реализацию технологии замкнутого топливного цикла.

Об авторах

А. А. Горюнов

Уральский федеральный университет имени первого Президента России Б.Н. Ельцина

Е. А. Тихонов

Уральский федеральный университет имени первого Президента России Б.Н. Ельцина

А. А. Власов

Уральский федеральный университет имени первого Президента России Б.Н. Ельцина

А. А. Бушуев

Уральский федеральный университет имени первого Президента России Б.Н. Ельцина

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

Статистика

Просмотры

Аннотация - 54

PDF (Russian) - 31

Ссылки

• Ссылки не определены.