CC BY
135
БЛОХИНА Александра Николаевна
инженер-конструктор отдела «Конструкторско-расчетный отдел прочности» (ОАО ОКБ «ГИДРОПРЕСС»)
СТОЛОТНЮК Сергей Владимирович
доцент, кандидат технических наук
СТОЛОТНЮК Ярослава Донатовна
доцент, кандидат технических наук
кафедры «Ядерные реакторы и установки» (МГТУ им. Н.Э. Баумана)
УДК 621.039.54
Сравнительный анализ напряженно-деформированного состояния конструкций тепловыделяющих сборок реактора ВВЭР-1000 при динамических нагрузках
А.Н. Блохина, С.В. Столотнюк, Я.Д. Столотнюк
Разработаны конечно-элементные модели тепловыделяющих сборок УТВС, ТВС-2 (конструкции ОКБ «ГИДРОПРЕСС») и ТВСА (конструкция «ОКБМ АФРИКАНТОВ» г. Нижний Новгород) в среде программного комплекса ANSYS. Определены собственные частоты и формы колебаний конструкций. Проведен расчет отклика конструкций на сейсмическое воздействие линейно-спектральным методом. Выполнен сравнительный анализ напряженно-деформированного состояния. Выявлены особенности предложенных конструкций.
Ключевые слова: тепловыделяющая сборка, конечно-элементная модель, напряженно-деформированное состояние, собственные частоты, формы колебаний, сейсмическое воздействие.
The finite-element models of the fuel assemblies UTVS, TVS-2 (designed by OKB «GIDROPRESS») and TVSA (designed by «AFRIKANTOV OKBM», N. Novgorod) were created by ANSYS code. The structure's eigene frequencies and oscillation modes were determined. The computation of the structure's response to seismic loads was based on the linear-spectral method. The comparative analysis of the stress-strained states was fulfilled. The features of the structure under consideration were discovered.
Keywords: fuel assembly, finite-element model, stress-strained state, eigene frequencies, oscillation modes, seismic loads.
Т_Геобходимость увеличения глубины выгорания ядерного топлива до величин более 55 МВт-сут/кг урана и перевода активных зон на четырех- и более годичные кампании [1] предполагает повышение надежности конструктивных элементов активных зон, главными из которых являются твэлы и органы регулирования. В настоящее время в нашей стране и за рубежом актуальной является проблема непроектного срабатывания аварийной защиты при увеличении глубины выгорания, заключающаяся, например, в превышении времени ввода регулирующих стержней (задержка в направляющих каналах), а в некоторых случаях недоход стержня до конечного положения. Эти случаи показали необходимость обоснования безопасности реакторной установки с точки зрения обеспечения эффективности аварийной защиты.
Подразумевается доведение режима аварийной защиты до проектных аварийных сценариев с максимальным расчетным землетрясением.
Представляется актуальной проблема сделать экономически приемлемыми компенсирующие меры проектирующих организаций, направленные на исключение застревания органов регулирования (ОР) системы управления и защиты (СУЗ). Минимизация искривления направляющих каналов ОР СУЗ привела к разработке усовершенствованных тепловыделяющих сборок (ТВС) более жесткой конструкции в ОКБ «ГИДРОПРЕСС» и «ОКБМ АФРИКАНТОВ» (г. Нижний Новгород). Необходимы эффективные методики оценки принимаемых конструктивных решений, учитывающие чувствительность конструкций ТВС к большим динамическим воздействиям.
В связи с увеличением сроков службы АЭС возрастает вероятность реализации таких внешних динамических воздействий как землетрясение, ударные волны техногенной природы, ураганы. Эти факторы существенно влияют на безопасность АЭС и требуют их учета многими нормативными документами. Сейсмические воздействия приводят к возникновению значительных инерционных нагрузок, отличающихся частотным составом и длительностью. Неожиданность реализации и значительная интенсивность таких динамических нагрузок предполагает всесторонний анализ прочности. Напряжения, обусловленные действием динамических нагрузок, часто являются определяющими при обосновании прочности.
Расчеты на динамические нагрузки — отдельная область прочностных расчетов, проводимых для обоснования конструкции с точки зрения безопасности. Расчетное моделирование ТВС может быть осуществлено на основе конечно-элементных технологий с учетом различных уровней приближений в описании геометрии элементов конструкции.
Данная работа посвящена разработке конечно-элементных моделей модернизированных ТВС реактора ВВЭР-1000 для сравнительного анализа влияния принимаемых конструктивных решений на повышение устойчивости
к формоизменению в условиях действия экстремальных нагрузок.
Разработка конечно-элементных моделей ТВС
В ходе выполнения работы были созданы конечно-элементные динамические модели усовершенствованной ТВС (УТВС) и ТВС-2, разработанной в ОКБ «ГИДРОПРЕСС», а также ТВСА, разработанных в «ОКБМ АФРИКАНТОВ» (рис. 1) и проведены следующие расчеты: собственных частот и форм колебаний кассет (модальный анализ); отклика на возмущение, характерное максимальному проектному землетрясению (спектральный анализ).
Y
б
□ || И I I I I II П I II I ¡Ш
я
Рис. 1. Конечно-элементные модели: а — УТВС; б — ТВС-2; в — ТВСА
Современный программный конечно-элементный комплекс ANSYS [2] позволяет построить сложные трехмерные динамические модели, учитывающие особенности реальной конструкции, повышающие точность расчетов и, в ряде случаев, снижающие излишний консерватизм. Однако описание локальных эффектов деформирования в сложной конструкции ТВС ведет к появлению значительного числа степеней свободы дискретной модели и невозможности получения численного результата. Поэтому для понижения числа степеней свободы в данной работе была использована упрощенная гипотеза о представлении дис-танционирующих и опорных решеток в виде упругих и жестких связей.
Существенным упрощением в этих расчетах является не учет влияния теплоносителя, как среды, демпфирующей колебания, а также не учет влияния подкрепляющего действия со
стороны соседних ТВС в составе активной зоны реактора. Данная модель не учитывает изменение механических свойств материалов кассеты и неравномерность их распределения вследствие воздействия тепловых потоков и потоков ионизирующих излучений.
Все представленные ТВС имеют схожее конструктивное исполнение [1, 3] и состоят из 18 направляющих каналов, 312 твэлов, центральной трубки, 15 дистанционирующих и одной опорной решеток, головки и хвостовика. Особенностью конструкции рассматриваемого пучка твэлов УТВС является наличие ограничительных втулок на трех направляющих каналах, которые в процессе эксплуатации жестко соединяются с дистанционирующими решетками (ДР). Это состояние можно считать потенциально опасным. В созданной конечно-элементной модели УТВС подобные соединения учитываются как сварные. ТВС-2 характеризуется наличием жесткого каркаса, сформированного направляющими каналами (НК), привариваемыми к дистанционирую-щим решеткам. В ТВСА каркас жесткости состоит из НК и уголков, приваренных к ДР.
Также при создании конечно-элементных моделей учитывались и другие особенности модернизированных ТВС: увеличенная высота и толщина ячеек ДР; усиленная головка и крепление хвостовика к опорной решетке в ТВС-2; в ТВСА усиленные НК между пучком твэлов и головкой.
Свойства материала головки и хвостовика соответствуют стали 08Х18Н10Т, остальных элементов — цирконию при рабочих температурах, характерных для реактора ВВЭР-1000 (325 °С). В качестве свойств наполняющего твэлы материала была задана плотность диоксида урана при тех же температурах.
Элементы ТВС цилиндрической формы (твэлы, направляющие каналы, центральная трубка, обечайка головки и хвостовик) моделируются с помощью конечного элемента pipe16 — упругая прямая труба. Данный элемент поддерживает растяжение, сжатие, изгиб и кручение, при этом задается двумя узлами и имеет шесть степеней свободы в каждом узле. Исходные данные элемента также включают
внешний диаметр трубы, толщину ее стенки и свойства наполняющего материала.
Дистанционирующие решетки, опорные решетки ТВС и решетки головки моделируются балочным элементом Ъеаш4 — трехмерная упругая балка. Элемент определяется двумя узлами, каждый из которых имеет по шесть степеней свободы, площадью поперечного сечения, двумя моментами инерции сечения и двумя толщинами. Элемент может воспринимать растяжение, сжатие, кручение и изгиб. Элементы жестко соединяют 331 узел решетки. При этом решетки представляют собой абсолютно жесткие конструкции, модуль упругости для которых был задан большим на четыре порядка по сравнению с соответствующими характеристиками для основных материалов (стали и циркония).
Для моделирования каркасных уголков ТВСА использовались элементы 8Ье1143 — оболочка со свойствами пластичности. Данный элемент имеет шесть степеней свободы в каждом узле. Задается четырьмя узлами и толщиной.
Для определения коэффициента жесткости в соединении твэл — ДР были отдельно построены оболочечные упрощенные конечно-элементные модели фрагментов ДР двух типов: для УТВС и модернизированных УТВС (рис. 2). Модели состоят из 19 ячеек, соединенных в форме шестигранника. Ячейки ДР моделировались конечными элементами типа 8Ье1143. Материал ячеек — цирконий. Соединение ячеек между собой моделирует контактную сварку по двум точкам и выполнено посредством соединения соответствующих соседних узлов ячеек абсолютно жесткими связями. Система закреплена по верхнему ряду ячеек.
Определение жесткости проводилось для центральной ячейки из заданного набора. К линии контакта с твэлом прикладывалась нагрузка, проводился расчет соответствующих перемещений, жесткость ячейки в составе решетки определялась, как сила, отнесенная к перемещению и оказалась равна 322,6 кН/м для УТВС и 770 кН/м для модернизированных УТВС.
Связь отдельных элементов конструкции между собой осуществлялась при помощи эле-
Рис. 2. Оболочечная модель дистанционирующей решетки
мента ша11г1х27 — матрица жесткости. Элемент не имеет определенной формы, но его упругий кинематический отклик может быть определен коэффициентами жесткости. Таким образом, элемент задается двумя узлами и коэффициентами матрицы жесткости. Жесткость подбиралась в соответствии с условиями скрепления элементов в рассматриваемых ТВС.
Коэффициенты жесткости задавались по трем направлениям линейных перемещений и трем направлениям угловых перемещений для пары узлов элемента ша11г1х27, перемещающихся под действием соответствующих сил и моментов.
Всего в моделях ТВС было использовано три типа связей: абсолютно жесткая связь, связь между узлами твэлов (и НК для УТВС и ТВС-2) и ДР, связь, моделирующая пружинный блок на головке ТВС. Для первого типа связи (абсолютно жесткая связь, моделирующая сварные и резьбовые соединения, а также цанговое крепление твэлов в опорной решетке) все коэффициенты жесткости, препятствующие смещениям узлов по направлениям шести степеней свободы, были заданы равными (по абсолютной величине) 1-108 Н/м.
Для второго типа связей (соединение твэл — ДР, НК — ДР) в поперечном направлении ТВС
были использованы коэффициенты жесткости, определенные из расчета отдельных моделей дистанционирующей решетки (322,6 кН/м для УТВС и 770 кН/м для ТВС-2 и ТВСА).
Для связи, моделирующей пружинный блок на головке ТВС, был задан коэффициент жесткости в осевом направлении для каждой из 19 пружин, равный 21 052,6 Н/м для УТВС и 17 710,5 Н/м для ТВС-2 и ТВСА. В поперечном направлении жесткость каждой пружины определялась жесткостью направляющего канала в этом направлении и составила 8 264,5 Н/м для всех моделей ТВС.
Модели были закреплены в узлах верхней обечайки головки и нижнего конца хвостовика ТВС. В верхней части кассеты было реализовано закрепление типа шарнир, при этом у головки остается возможность осевого перемещения. Хвостовик закреплен абсолютно жестко (реализовано закрепление типа заделка).
Верификация моделей
По результатам статических испытаний и расчетов, проведенных в ОКБ «ГИДРОПРЕСС» [4], были верифицированы модели УТВС и ТВС-2. При этом зависимость перемещений УТВС и ТВС-2 от приложенной осевой силы в диапазоне до 120 кН практически совпали. При нагрузке поперечной силой зависимость максимального прогиба ТВС-2 совпадает с расчетными и экспериментальными зависимостями ОКБ «Гидропресс» с точностью до погрешности.
Также модель ТВС-2 была верифицирована по результатам расчета собственных частот, определенных по коду И/ОЯ в РНЦ «Курчатовский институт» («КИ») [5]. Сравнение результатов расчетов приведено в табл. 1.
Таблица 1
Расчетные значения собственных частот колебаний
Расчетная модель Собственные частоты изгибных форм, Гц
Созданная модель ТВС-2 3,61 7,71 12,66 17,94 23,77
По коду и/ОЯ РНЦ «КИ» 3,38 7,51 11,64 17,025 24,05
Результаты проведенных расчетов показывают хорошее совпадение с экспериментальными
и полученными с использованием кода Ц^ОЯ данными. Различие состоит в пределах погрешности экспериментальных и расчетных данных.
Использование конечных элементов, близких по форме, позволило получить в моделях регулярную сетку разбиения. С такой сеткой можно получить сравнительно точные результаты при использовании конечных элементов с малым числом узлов, что позволяет существенно сократить объемы вычислений и сэкономить машинное время.
На практике адекватность конечно-элементной сетки оценивается значениями нерав-номерностей распределения расчетной величины (напряжений, смещений, температур) по соседним элементам. В программе ЛМ8У8 используются такие оценки как: относительная погрешность между результатами в точках интегрирования элементов и результатами, усредненными по узлам; оценки отклонения вероятной нижней и верхней границ результатов относительно вычисленных минимальных и максимальных величин.
Для исследуемых в данной работе моделей ТВС, которые выполнены с использованием регулярной сетки, значения перечисленных погрешностей варьируются в пределах 3...5%, что позволяет сделать вывод об адекватности полученных с использованием этих моделей результатов.
Расчет отклика ТВС на сейсмическое воздействие
Расчет отклика ТВС на сейсмическое воздействие в поперечном направлении требует предварительного определения собственных частот и форм колебаний кассеты в соответствующем направлении. Указанные параметры колебаний кассеты определяются из модального анализа, в ходе которого решаются системы линейных уравнений:
[М]{ы"} + [ К\{ы} = {0}. (1)
Матричное уравнение (1) является частным случаем общего уравнения движения произвольной системы:
[М ]{ы"} + [С ]{ы'} + [ К\{и} = {^ ()}, (2)
где [М] — матрица масс системы; [С] — матрица сопротивлений; [К — матрица жесткости; {ДО} — функция нагрузки, зависящая от времени; {и} — вектор узловых перемещений, {и'} — вектор узловых скоростей; {и''} — вектор узловых ускорений.
При определении собственных частот и форм колебаний конструкции предполагается, что совершаются свободные незатухающие колебания, т. е.
{^(Г)} = {0} и [С\ = [0\. (3)
В спектральном анализе искомые узловые перемещения находят из решения уравнения движения системы в общем виде.
Для перехода от узловых смещений к величинам напряжений используется зависимость между перемещениями и деформациями, а также зависимость между деформациями и напряжениями (закон Гука):
{в} = [Д\{и}; {о} = [Е\{в}, (4)
где [Д| — матрица операторов дифференцирования; [Е] — матрица, членами которой являются модуль упругости и модуль сдвига.
В соответствии с рекомендациями, приведенными в «Нормах расчета на прочность...» [6], расчет ТВС на сейсмическое воздействие был проведен линейно-спектральным методом. При данном методе расчета нагрузки, имитирующие максимальное проектное землетрясение, задаются в виде обобщенных спектров ответа горизонтальных колебаний, другими словами, к опорам конструкции прикладывается спектр ускорений в зависимости от частоты колебаний. Вид спектра ускорений и численные значения взяты из «Норм расчета на прочность.» [6] и соответствуют сейсмичности площадки 9 баллов по шкале М8К-64 для конструкций, располагаемых на уровне земли, при относительном демпфировании колебаний 0,02.
Нагрузка, действующая в узлах конструкции, при линейно-спектральном методе суммируется по всем определенным ранее формам колебаний по формуле
а и в
Рис. 3. Изгибные формы свободных незатухающих колебаний системы: а — УТВС; б — ТВС-2; в — ТВСА
NK =
2 Ч
j =1
(5)
где Nк — расчетное усилие в К-м узле системы; S — номер максимальной собственной частоты (формы), определенной при модальном анализе; у — номер формы колебаний.
Результаты расчетов
В результате модального анализа были получены изгибные и крутильные формы свободных незатухающих колебаний системы с соответствующими значениями собственных частот (рис. 3, 4).
Результаты расчета ТВС на сейсмическое воздействие (табл. 2) представлены в виде мак-
симальных смещений конструкции от положения равновесия в виде формы колебаний и значений напряжений, действующих в каркасе жесткости конструкции.
Наибольшие напряжения в каркасе жесткости УТВС возникают в трех направляющих каналах, приваренных к дистанционирующим решеткам. При этом максимум напряжений приходится на участок НК между твэльным пучком и головкой.
Поскольку в ТВС-2 все НК приварены к ДР, распределение напряжений в них одинаково и меньше, чем в подобных узлах УТВС. Максимальное значение также приходится на участок НК между твэлами и головкой, но абсолютное значение эквивалентных напряжений значительно меньше чем в УТВС.
Рис. 4. Крутильные формы свободных незатухающих колебаний системы: а — УТВС; б — ТВС-2; в — ТВСА
Таблица 2
Результаты расчетов
Величина УТВС ТВС-2 ТВСА
Собственные частоты изгиб-ных колебаний по направлению оси Оу (грань шестигранника), Гц 2,12 4,61 7,73 11,63 16,49 22,44 3,61 7,71 12,66 17,94 23,77 3,27 7,47 12,52 17,96 24,09
Собственные частоты изгиб-ных колебаний по направлению оси Оx (угол шестигранника), Гц 2,14 4,68 7,85 11,82 16,76 22,8 3,64 7,8 12,81 18,15 24,05 3,32 7,58 12,68 18,21 24,41
Собственные частоты крутильных колебаний, Гц 4,99 12,14 19,81 8,15 16,33 24,58 24,81
Окончание табл. 2
Величина УТВС ТВС-2 ТВСА
Максимальный прогиб, мм 50,59 19,97 23,19
Максимальные эквивалентные напряжения в НК, МПа 371 168 82,4
Максимальные напряжения в каркасных уголках, МПа — — 226
Применение уголков в каркасе жесткости ТВСА, а также специальных стальных втулок, приваренных к НК в наиболее ослабленном месте каркаса (между твэлами и головкой), позволило значительно разгрузить направляющие каналы, являющиеся основной частью каркаса жесткости. Максимальные напряжения в НК снизились по сравнению с ТВС-2 примерно в 2 раза.
Максимальные эквивалентные напряжения в каркасах жесткости представлены на рис. 5—7.
Рис. 5. Максимальные эквивалентные напряжения в каркасе жесткости УТВС
Рис. 6. Максимальные эквивалентные напряжения в каркасе жесткости ТВС-2
Рис. 7. Максимальные эквивалентные
Сравнительный анализ напряженно-деформированного состояния различных типов ТВС показал, что модернизированные ТВС-2 и ТВСА позволяют значительно снизить напряжения в направляющих каналах по сравнению с УТВС. Направляющие каналы, используемые в рассмотренных ТВС для размещения ОР СУЗ и образующие каркас жесткости конструкций, определяют эффект формоизменения, который за время кампании не должен приводить к нарушениям в работе ОР. Модернизированные конструкции ТВС позволяют перенести часть напряженного состояния к дистационирующим решеткам и каркасным уголкам. Если рассматривать ТВСА в целом, наиболее напряженные участки конструкции расположены в местах сварного соединения каркасных уголков с ободами дистанциони-рующих решеток (рис. 8). Увеличение нагрузки в новых конструкциях дистанционирующих решеток модернизированных ТВС позволяет значительно разгрузить каркас жесткости.
напряжения в каркасе жесткости ТВСА
В ТВС-2 напряжения в направляющих каналах уменьшаются в 2 раза, а в ТВСА — в 4 раза.
Расчет напряженно-деформированного состояния в ячейках дистанционирующих решеток модернизированных ТВС показал, что напряжения возрастают от центральных ячеек к периферийным.
В ячейках дистанционирующих решеток ТВС-2 и УТВС максимальные напряжения возникают в ячейках с твэлами, прилегающих к середине грани сварного бандажа дистанцио-нирующей решетки. В ТВСА максимальные напряжения возникают в угловых ячейках, в местах сварного соединения каркасных уголков с ободами дистанционирующих решеток (рис. 9). Учитывая, что в межкассетном зазоре возле этих ячеек генерируется максимум пуль-сационной энергии потока возможно ожидать в указанных ячейках повышенный фреттинг износ оболочек твэлов.
Результаты расчетов показали, что повышение жесткости каркаса УТВС привело к увеличению
Рис. 9. Распределение эквивалентных напряжений в дистанционирующих решетках ТВС:
а — УТВС; б — ТВСА
значений собственных частот колебаний. Собственные частоты изгибных форм колебаний модернизированных ТВС увеличились примерно в 1,5 раза. При этом количество первых изгибных форм колебаний, соответствующих реализации пиковых значений ускорений при сейсмическом воздействии в низкочастотной области (до 15 Гц), уменьшилось с 4 до 3. Частоты крутильных форм колебаний в ТВС-2 значительно увеличиваются, а в ТВСА частота первой формы оказалась равной 24 Гц, что более чем в 4 раза выше, чем у УТВС. Применение уголков в каркасе жесткости ТВСА привело к значительному увеличению низшей собственной частоты крутильных колебаний. Повышение частоты колебаний и уменьшение числа форм должно благоприятно сказаться на условиях нагружения и износа поверхности оболочек твэла при взаимном угловом перемещении ДР-твэл.
Созданием новых конструкций ТВС усилено подкрепляющее действие всего массива ТВС активной зоны. Значительно снижен возможный максимальный прогиб ТВС, что на ряду с уменьшением зазоров между ТВС должно привести к меньшему изменению продольных осей направляющих каналов при сейсмических нагрузках.
Расчетный анализ отклика конструкции ТВС на сейсмические нагрузки показал, что изменение конструкции ТВС от УТВС к ТВС-2 и ТВСА привело к снижению напряжений в элементах каркаса жесткости. Повышенная жесткость конструкции ТВС-2 и ТВСА обеспечивает вдвое меньшие максимальные прогибы при заданном возмущении, чем у базовой УТВС.
В ТВС-2 приварка дистанционирующих решеток к направляющим каналам делает конструкцию более устойчивой к прогибам, чем конструкция с приваркой дистанционирующих решеток к уголкам в ТВСА.
В ТВСА в местах расположения направляющих каналов напряжения в ячейках ДР значительно ниже чем в периферийных, что должно позволить снизить трудности с движением по-
глощающих стержней (ПС) СУЗ во время динамических нагрузок. На условия нагружения направляющих каналов в ячейках ДР оказывает также влияние продольное поджатие пружинами головки ТВС. Требуется анализ совместного нагружения направляющих каналов осевыми и динамическими нагрузками. Дальнейшее расчетно-экспериментальное обоснование предложенных конструкций ТВС с учетом сопряженных нагрузок позволит проанализировать влияние изменений условий поджатия ТВС в течении кампании на напряженно-деформированное состояние элементов конструкций.
Методика-технология проведения расчетного анализа реализована в виде программ, позволяющих использовать ее в инженерной практике на стадиях разработки новых модификаций ТВС.
Использованная методика-технология создания расчетных динамических моделей сложных конструкций на примере ТВС реактора ВВЭР-1000 позволяет проводить прочностное обоснование работоспособности при экстремальных динамических нагрузках [7].
Литература
1. Реакторы ВВЭР-1000 для атомных электростанций. В.К. Резепов, В.П. Денисов, Н.А. Кирилюк и др. М.: ИКЦ «Академкнига», 2004. 333 с.
2. Каплун А.Б. ANSYS в руках инженера. М.: Едиториал УРСС, 2004. 272 с.
3. Шмелев В.Д., Драгунов Ю.Г., Денисов В.П., Васильчен-ко И.Н. Активные зоны ВВЭР для атомных электростанций. М.: ИКЦ «Академкнига», 2004. 220 с.
4. Шарый Н.В., Семишкин В.П., Пиминов В.А, Драгунов Ю.Г. Прочность основного оборудования и трубопроводов реакторных установок ВВЭР. М.: ИздАТ, 2004, 496 с.
5. Киселев А.С., Киселев А.С., Даничев В.В. Аннотация программы «UZOR» // Вант. Сер. Физика ядерных реакторов. 1999. Вып. 1. С. 109—113.
6. Нормы расчета на прочность оборудования и трубопроводов атомных энергетических установок. ПНАЭГ-7-002—86. М.: Энергоатомиздат, 1989. 525 с.
7. Прочность и надежность конструкций АЭС при экстремальных воздействиях/ С.Е. Бугаенко, С.Л. Бугорин, ГС. Шуль-ман, С.Г. Шульман. М.: Энергоатомиздат, 2005. 576 с.
Статья поступила в редакцию 22.03.2012
