CC BY
43
СОРОКИН Федор Дмитриевич
профессор
ЕСИН Максим Юрьевич
аспирант кафедры «Прикладная механика»
ПЕРЕВЕЗЕНЦЕВ Владимир Васильевич
кандидат технических наук, доцент кафедры «Ядерные реакторы и установки» (МГТУ им. Н.Э. Баумана) e-mail: [email protected]
УДК 621.039.524: 532.5, 534
Моделирование связанных гидроупругих колебаний тепловыделяющих сборок в активной зоне ВВЭР-440
Ф.Д. Сорокин, М.Ю. Есин, В.В. Перевезенцев
Для моделирования связанных гидроупругих колебаний тепловыделяющих сборок (ТВС) 7-кассетного фрагмента активной зоны реактора ВВЭР-440 вычислены гидродинамические коэффициентов влияния ТВС друг на друга. Задача гидродинамики решалась в плоской постановке методом конечных элементов. Определение равнодействующих сил давления воды на ТВС при гармонических колебаниях одной из кассет позволило последовательно определить гидродинамические коэффициенты влияния (матрица 14x14) и полную матрицу масс. В результате исходная задача сведена к традиционной задаче динамики многомассовых систем.
Ключевые слова: ядерные реакторы ВВЭР, гидроупругие колебания ТВС, гидродинамические коэффициенты влияния, частоты и формы собственных колебаний.
Modeling of bound hydrodynamic vibrations of fuel assemblies in WWER-440 reactor core
F.D. Sorokin, M.Yu. Esin, V.V. Perevezentsev
The hydrodynamic coefficients of fuel assemblies influence on each other have been calculated for modeling the bound hydrodynamic vibrations in a fragment of WWER-440 reactor core composed of seven fuel assemblies. The hydrodynamic problem has been solved as a plane problem of FEM. Calculating of the resultant of water pressure on the fuel assemblies during the harmonic motion of one fuel assembly allowed to sequentially determine the hydrodynamic coefficients of influence (14x14 matrix) and the total matrix mass. As a result, the initial problem is reduced to the traditional problem of the dynamics of multimass systems.
Keywords: water-water power reactor (WWER), hydrodynamic vibrations of fuel assemblies, hydrodynamic coefficients of influence, frequency and form of natural vibrations.
ТЭибрационные процессы в ядерных реакторах в значительной степени определяют надежность оборудования и безопасность эксплуатации энергоблока в целом. Динамическое воздействие потока теплоносителя на конструктивные элементы реакторных систем инициирует и поддерживает их механические колебания (вибрации). Одними из наиболее важных механических систем, определяющих
безаварийную работу реакторных установок ВВЭР, являются тепловыделяющие сборки (ТВС). При эксплуатации реакторных установок ВВЭР наблюдались повреждения ТВС, обусловленные вибрационными процессами [1, 2]. Выявление механизмов возникновения вибраций ТВС в реакторных установках ВВЭР — важнейший этап создания методик расчетных и экспериментальных анализов вибрационных процессов и выработки научно обоснованных рекомендаций по конструктивным решениям, направленным на повышение устойчивости ТВС в целом и ее отдельных элементов к динамическим воздействиям теплоносителя. Для разработки расчетных моделей связанных гидроупругих колебаний ТВС в активной зоне реактора необходимы данные по динамическим характеристикам (собственным частотам и формам колебаний) гидроупругой системы.
В данной работе предлагается методика численного расчета собственных частот и форм связанных гидроупругих колебаний кассет 7-кассетного фрагмента активной зоны ВВЭР-440 (рис. 1).
В общем виде гидроупругие свободные колебания системы одинаковых тел, погруженных в жидкость, описываются следующая системой уравнений:
Рис. 1. Поперечное сечение 7-кассетного фрагмента активной зоны ВВЭР-440:
1 — ТВС; 2 — заполненный водой зазор между ТВС
([M] + [ AM]){q} + [C]{q} = 0;
"1 0 0 0' 0 10 0
[M] =
m,
[C] = .
0 0 ... 1
10 0 0
0 10 0
0 0 ... 1
(1)
где {ц} — вектор узловых перемещений средних сечений кассет; [М] — матрица приведенных масс без учета жидкости; [ДМ] — матрица гидродинамических коэффициентов влияния; [С] — матрица приведенных жесткостей. Крутильные степени свободы не учитывались, поэтому число степеней свободы равно 14 (по два линейных перемещения у каждого из семи средних сечений ТВС).
По экспериментальным данным о собственных частотах изгибных колебаний кассет в воздухе удалось оценить изгибную жесткость конструкции. Поэтому для расчета собственных частот и форм системы ТВС, погруженных в объем воды в активной зоне достаточно вычислить присоединенные массы и гидродинамические коэффициенты влияния ТВС друг на друга, т. е. решить задачу гидродинамики.
Для решения задачи гидродинамики в постановке Лагранжа применялся конечный элемент жидкости FLUID79 из библиотеки комплекса ANSYS, который является модификацией 4-узлового элемента плоской задачи теории упругости. Элемент предназначен для решения плоских задач гидростатики и взаимодействия идеальной сжимаемой жидкости с твердыми телами. Аналогичный элемент был также реализован в специализированном конечно-элементном комплексе на основе методик, изложенных в работе [5].
Функции формы конечного элемента определяют обычными билинейными выражениями, которые используют при решении плоской задачи теории упругости. Соотношения между деформациями и напряжениями в жидкости имеют вид
I ху Юz
К 0 »
» I « 0 0 В
р
м.
(2)
где вь — объемная деформация; уху — угол сдвига; юz — малый поворот; р — гидростатическое давление; тху — касательное напряжение; Мг — фиктивный скручивающий момент.
Деформации и повороты связаны с перемещениями обычными соотношениям механики сплошной среды [6]:
дЫ дУ
в ь = в х + в у =--\--;
ь х у дх ду
_ дЫ дУ
Уху ~~ду +
Ю,
Величины тху и Мг необходимо вводить в рассмотрение только для того, чтобы элемент обладал некоторой сдвиговой и вращательной устойчивостью, т. е. чтобы не возникало ошибок при численном счете.
Единственной существенной физической постоянной (помимо плотности жидкости) для данного элемента является модуль объемной упругости жидкости К (для воды он составляет порядка 2-109 Па). Остальные константы, входящие в уравнение состояния, задаются на 9 порядков меньшим числом, чем К:
£ = В = 10-9 К. (4)
Фактически из (2)—(4) с учетом большой величины модуля объемной упругости следует традиционное уравнение состояния неизменности объема несжимаемой жидкости:
дЫ дУ „ — \ — = 0.
дх ду
Константы S и В вводятся по той причине, что из соображения численного характера нельзя полностью отказаться от восприятия конечным элементом сдвигов и вращений, не-
смотря на то, что в модели идеальной сжимаемой жидкости таких констант не существует. Подчеркнем, что константа S в формуле (2) имеет смысл фиктивной сдвиговой жесткости и поэтому не имеет никакого отношения к вязкости (жидкость идеальная).
Матрица жесткости конечного элемента строилась энергетическим методом. При этом удельная энергия деформаций определялась выражением
=2 к [в 2 +(у 2у + ю| )кг' ].
и0 = 2
Очевидно, что в выражении (4) существенную роль играет только объемная деформация.
Матрица масс элемента вычислялась обычным образом через функции формы
[м(*) ] = рЯМ т №А,
где р — плотность жидкости; [ЭД — матрица функций формы; А — площадь элемента.
Следует отметить, что в элементе FLUID79 из ANSYS применяется диагональная матрица масс, в авторской реализации конечного элемента матрица масс была полностью заполненной.
Конечный элемент тестировался на нескольких задачах, одной из которых была классическая задача о колебаниях коаксиально расположенного цилиндра в жидкости (задача Стокса) (рис. 2), имеющая аналитическое решение [3].
В задаче Стокса длинный жесткий цилиндр совершает вынужденные гармонические колебания с частотой / в длинной цилиндрической полости, заполненной жидкостью. Были приняты следующие числовые значения для параметров данной системы: а = 0,2 м, кг
Ь = 0,3 м, р = 1 000 —г, / = 3 Гц. Полость была м3
разбита на 1 200 конечных элементов (по радиусу — 10 элементов, по окружности — 120 элементов). На границе неподвижного цилиндра перемещения запрещались. На границе подвижного цилиндра были заданы перемещения, связанные с амплитудой центра цилиндра:
в
ь
т
ху
А
Рис. 2. Тестовая задача о колебаниях коаксиально расположенного цилиндра в жидкости (задача Стокса)
u = 0;
v = A ео8 ,
(7)
где u — перемещение по оси х0; v — перемещения по оси у0; ш = 2п/.
То, что перемещения меняются по гармоническому закону (7), учитывалось опциями комплекса ANSYS. Присоединенная масса на единицу длины цилиндра вычислялась через равнодействующую сил давления по формуле
Дm
1 2 п
= " АШ2 ( P о(ф)'
й ф^ф,
(8)
где р0 — амплитудное значение давления.
Результаты определения присоединенной массы для тестовой задачи с использованием комплекса ANSYS (МКЭ), аналитическим [3] и численно-аналитическим [4] методами оказались достаточно близкими по значению (табл. 1).
Таблица 1
Сопоставление результатов расчета различными методами присоединенной массы для тестовой задачи колебаний коаксиального цилиндра
Метод определения Аналитический [3] Числен-но-аналити-ческий [4] МКЭ
Удельная присоединенная масса Дт, кг/м 383,6 384 383,7
При использовании численно-аналитического метода [4] необходимо учитывать вязкость воды, для которой принималось обычное значение ц = 0,001 Па с. Как видно из табл. 1, влияние вязкости жидкости в рассматриваемой задаче практически не заметно. Таким же незначительным оказалось влияние частоты на величину Дт, если даже частота меняется на порядок.
Колебания в жидкости коаксиального цилиндра (задача Стокса) существенно отличаются по геометрии от 7-кассетного фрагмента активной зоны ВВЭР-440. В связи с этим более корректное тестирование конечного элемента выполнялось на примере 3-кассетного фрагмента активной зоны (рис. 3), для которого может быть применен численно-аналитический метод [4].
Рис. 3. Пример 3-кассетного фрагмента активной зоны, колеблется центральная кассета
При исследовании 3-кассетного фрагмента активной зоны были приняты те же предположения, что и в задаче Стокса для коаксиального цилиндра. Числовые значения параметров 3-кассетного фрагмента были заданы соответствующими реальным размерам кассет и плотности воды в активной зоне ВВЭР: а =
кг
= 0,08 м, к = 0,0055 м, р = 1 000 —; значение
м3
эффективного зазора к выбиралось в соответствии с работой [4].
Для применения численно-аналитического метода согласно [4] достаточно найти распределение средней скорости жидкости по границам ТВС. Как видно на рис. 3 существенным отличием задачи 3-кассетного фрагмента от задачи Стокса является возможность разветвления потока в узлах А, В, С, Б. Естественно
предположить, что поток в узле А (и соответственно в узле С) делится в соответствии с гидравлическими сопротивлениями путей, по которым течет жидкость. Путь А—В в 5 раз меньше, чем обходной путь вокруг бокового шестиугольника А—Ь—В, поэтому поток в направлении А—В должен быть в 5 раз больше, чем поток по обходному пути А—Ь—В.
Амплитудное значение средней скорости жидкости, втекающей в точку А со стороны ребра ОА равно отношению потока, вытесняемого гранью шестиугольника ОА, к ширине зазора
Аюа . п
^т- (9)
У ср1 =-
к
3
Тогда, в соответствии со сказанным выше, вызванные колебаниями ТВС скорости течения жидкости определяются следующими выражениями:
5
V ч = — V , =
ср2 6 ср1
1
V о = — V ,=■
ср3 6 р1
5Аюа . п
81П3
Аюа . п
6к
(10)
1-В; уср3 —
где Уср2 — скорость в направлении А— скорость в направлении А—Ь—В.
Соотношения (9), (10) записаны с учетом того, что поток в точке А разделяется в пропорции 5:1. Скорость жидкости в зазоре распределена симметрично относительно оси у0 (см. рис. 3). Найденное распределение амплитуды средней скорости по периметру канала показано на рис. 4.
Дальнейший расчет сил, действующих на шестиугольные сечения со стороны жидкости, полностью идентичен методике, приведенной в работе [4], с учетом другого закона распреде-
Рис. 4. Распределение амплитуды средней скорости потока:
а — по контуру О—А—В—М; б — по контуру А—Ь—В—А
ления скорости жидкости по контуру шестиугольника. В результате были найдены равнодействующие действительной и мнимой составляющих сил давления на центральную и боковые кассеты. Действительная составляющая сил давления, приложенная к центральной (подвижной) кассете, преобразуется в присоединенную массу аналогично (8). Действительные составляющие сил давления, действующие на боковые (неподвижные) кассеты преобразуется в гидродинамические коэффициенты влияния.
Та же задача о 3-кассетном фрагменте активной зоны решалась МКЭ. Процедура построения модели и задание граничных условий аналогична задаче Стокса. Конечно-элементная сетка выбиралась достаточно густой — отрезок по ширине зазора был разбит на 30 конечных элементов.
Конструкция нагружалась кинематически, для этого задавалось движение по гармоническому закону центрального узла средней кассеты. Остальные кассеты были закреплены неподвижно.
Наибольший интерес представляла проверка предположения о разделении потока в пропорции 1:5 в узле А. Амплитудные значения перемещений в области, примыкающей к этому узлу показаны на рис. 5. Как видно на рисунке, поток после прохождения узла А действительно делится в соотношении 1:5. Этим подтверждается предположение о разделении потока, принятое выше в численно-аналитическом методе.
Полученные результаты приведены в табл. 2. Результаты решения задачи двумя разными методами достаточно близки, что подтверждает надежность методики.
Таблица 2
Сравнение результатов расчета присоединенных масс и коэффициентов влияния численно-аналитическим методом [4] и МКЭ
ТВС Определяемые величины Численно-аналитический подход [4] МКЭ
Центральная Присоединенная масса Лотп, кг/м 172 168
Боковая Коэффициент влияния Лт12, кг/м 22,05 21
Рис. 5. Амплитуды перемещений жидкости в окрестности узла А
Таким образом, методика вычисления инерционных гидродинамических коэффициентов влияния для связанных колебаний ТВС в условиях 3-кассетного фрагмента активной зоны хорошо согласуется с численно-аналитическим методом [4].
Разработанная методика применялась для решения основной задачи о гидроупругих связанных колебаниях ТВС в 7-кассетном фрагменте активной зоны. Отличие от 3-кассетного фрагмента заключалось всего лишь в большем количестве степеней свободы. Среднее сечение каждой из семи кассет может перемещаться в двух направлениях (вращение кассет не учитывается), т. е. для всего фрагмента существует 14 степеней свободы. Каждая из указанных 14 степеней свободы по очереди приводилась в колебательное движение по гармоническому закону Давление, действующее по пе-
риметру сечений кассет, преобразовывалось в гидродинамические коэффициенты влияния аналогично (8) и [4]. Каждый из 14 запусков программы с одной активной степенью свободы дает 14 амплитуд сил, действующих на кассеты, т. е. позволяет заполнить одну строчку матрицы коэффициентов влияния [ДM]. После 14 запусков программы вся матрица коэффициентов влияния 14x14 оказалась заполненной.
Найденная матрица гидродинамических коэффициентов использовалась для расчета час-
тот и форм гидродинамически связанных колебаний кассет фрагмента.
При отсутствии жидкости [ДM] = [0] и, следовательно, собственные частоты всех кассет одинаковы. Удельная масса каждой кассеты и частота собственных колебаний в воздухе известны из наблюдений. Это позволяет найти приведенную (к среднему сечению) удельную жесткость кассеты
с = (2л/в)2 то = 7,74 • 104 Н/ м, (10)
где /в = 5,0 Гц — низшая частота собственных поперечных колебаний кассеты в воздухе; т0 = = 78,4 кг/м — погонная масса кассеты.
Таким образом, матрица масс ([M] + [ДM]) и жесткостей [^ в уравнении движения (1) определены. По этим значениям стандартными методами теории колебаний были найдены 14 собственных частот и собственных форм колебаний. Собственные частоты колебаний 7-кассетного фрагмента активной зоны представлены в табл. 3.
Таблица 3
Собственных частоты колебаний семикассетного фрагмента активной зоны
Номер формы Частоты собственных колебаний, Гц
1 1,964
2 1,971
3 1,984
4 2,122
5 2,134
6 2,159
7 2,166
8 2,168
9 2,748
10 2,748
11 2,877
12 3,897
13 3,899
14 4,241
Формы колебаний для 1- и 14-й частоты в виде векторные полей изображены на рис. 6.
ДЭ
щр
а б
Рис. 6. Собственные формы колебаний для двух частот:
а — / = 1,964 Гц; б — / = 4,241 Гц
Отметим, что в силу симметрии задачи собственные частоты, приведенные в табл. 3, располагаются группами из нескольких очень близких частот. Например, выделяется группа их трех первых частот, схожесть которых легко объяснить наличием трех равноправных направлений на рис. 6, а (очевидно, что поворот рис. 6, а на угол 60° приводит к другой форме, но с той же частотой). Исключение составляет последняя 14-я частота, не имеющая близких частот в таблице. Это объясняется тем, что поворот соответствующей ей формы (рис. 6, б) на 60°, приводит к той же самой форме (рисунок не меняется).
Сопоставление рассмотренных форм показывает, что в первом случае (см. рис. 6, а) движение жидкости стеснено намного сильнее, чем во втором (см. рис. 6, б), что приводит к гораздо большей присоединенной массе и, следовательно, к гораздо меньшей частоте. Для формы, представленной на рис. 6, б, жидкость как бы увлекается периферийными кассетами в круговое движение вокруг центральной кассеты, т. е. жидкость совершает движение близ-
кое к жесткому смещению. Этим объясняется малая присоединенная масса для 14-й формы колебаний.
Выводы
1. На основе МКЭ предложена методика расчета инерционных гидродинамических коэффициентов влияния для системы безчехло-вых кассет в макронеподвижной жидкости.
2. Методика оттестирована сопоставлением с численно-аналитическим методом [4].
3. По предложенной методике найдены гидродинамические коэффициенты влияния и вычислены частоты и формы собственных колебаний семикассетного фрагмента активной зоны.
Литература
1. Шмелев В.Д, Драгунов Ю.Г, Денисов В.П, Васильченко И.Н. Активные зоны ВВЭР для атомных электростанций. М.: ИКЦ «Академкнига», 2004. С. 220.
2. Денисов В.П., Драгунов Ю.Г. Реакторные установки ВВЭР для атомных электростанций. М.: ИздАТ, 2002. С. 387.
3. Синявский В.Ф., Федотовский В.С., Кухтин А.Б., Тере-ник Л.В. Инерционность и гидродинамическое демпфирование при колебаниях труб и трубных пучков жидкости // Динамические характеристики и колебания элементов энергетического оборудования. М.: Наука, 1980. С. 86—97.
4. Солонин В.И., Перевезенцев В.В., Сорокин Ф.Д. Демпфирование колебаний пучка твэлов чехловых тепловыделяющих сборок водоохлаждаемых реакторов // Вестник МГТУ им. Н.Э. Баумана. 2008. № 3. С. 75—85.
5. Zienkiewicz O.C., Taylor R.L. Finite Element Method. Vol. 1. The Basis. Amsterdam: Elsevier, 2000. P. 689.
6. Тимошенко С.П., Гудьер Дж. Теория упругости. М.: Наука, 1975. С. 560.
Статья поступила в редакцию 07.06.2012
