CC BY
121
ВШЗЕХЭШШ] выкэдшех ©аведжшй
СОРОКИН Федор Дмитриевич
профессор
КРУТЬКО Евгений Сергеевич
аспирант кафедры «Прикладная механика» (МГТУ им. Н.Э. Баумана) e-mail: sorokin fd@mail.ru
УДК 532.5, 534
Расчет присоединенной массы и коэффициента демпфирования для вибрирующего в цилиндрическом канале жесткого цилиндра на основе численного интегрирования уравнений движения вязкой жидкости
Ф.Д. Сорокин, Е.С.Крутько
Рассмотрено одномерное решение уравнений Навье — Стокса для расчета присоединенной массы и коэффициента демпфирования при вибрациях жесткого цилиндра в цилиндрическом канале. Гармонический характер колебаний и круговая геометрия канала позволяют исключить две независимые переменные, что приводит к обыкновенному дифференциальному уравнению 4-го порядка. Малая вязкость жидкости порождает узкие пограничные слои, вследствие чего полученное разрешающее уравнение является жестким. Указанная численная проблема преодолевается повышением точности (60—80 десятичных знаков в числах пакета Mathematica) и подбором шага интегрирования. Построены зависимости коэффициентов присоединенной массы и демпфирования от вязкости жидкости.
Ключевые слова: присоединенная масса, коэффициент демпфирования, цилиндрический канал.
Added mass and damping coefficient calculation for the rigid cylinder vibrating in cylindrical channel based on viscous fluid motion equation numerical integration
F.D. Sorokin, E.S. Krutko
The one-dimensional Navier-Stokes equations solution has been considered to perform the added mass and damping coefficient calculation during vibration of the rigid cylinder in the cylindrical channel. A harmonic type of oscillations and the cylindrical channel geometry allow to eliminate two independent variables, which leads to the fourth-order ordinary differential equation. Low liquid viscosity generates thin boundary layers, and consequently the derived resolving equation becomes stiff. This numerical problem can be solved by means of precision increasing (60—80 decimal digits in Mathematica software numbers representation) and appropriate integration step selection. The relations of added mass and damping coefficients on liquid viscosity have been plotted.
Keywords: added mass, damping coefficient, cylindrical channel.
Для расчета частот и амплитуд колебаний элементов энергетических установок в жидкости должны быть известны значения присоединенной массы и коэффициента демпфирования. Задача определения присоединенной массы и коэффициента демпфирования при колебаниях труб и трубных пучков в жидкости рассмотрена во многих экспериментальных и теоретических работах, в частности, в [1—3]. В работе [1] приведены приближенные аналитические выражения для расчета указанных величин, которые обобщены для трубных пучков различной конфигурации.
Целью данной статьи является расчет присоединенной массы и коэффициента демпфирования на основе точного (точнее, численно точного) решения уравнений Навье — Стокса. В отличие от приближенных методик предлагаемая методика позволяет исследовать случай взаимного влияния пограничных слоев, что имеет место при очень малом зазоре между цилиндрами, либо при очень высокой вязкости. Достоинством численных точных решений уравнений Навье — Стокса является то, что на их основе удается проводить тестирование и настройку современных конечно-элементных методик, например, правильно выбирать размер конечного элемента и т.п.
Рассмотрим движение жидкости в зазоре между двумя цилиндрами (рис. 1). Наружный цилиндр неподвижен, а внутренний совершает гармонические колебания по закону ^со8(ш?).
Уравнения движения вязкой несжимаемой жидкости для малых колебаний целесообразно рассматривать в постановке Лагранжа:
Эй
д2 и
д?
V-и = 0
= -Ур + цУ
д?
(1)
где и — вектор малых перемещений частиц жидкости относительно неподвижной системы отсчета; р — плотность; ц — динамическая вязкость; р — давление; ? — время. Постановка Лагранжа в рассматриваемой задаче более уместна, так как частицы жидкости не уходят в бесконечность, как в постановке Эйлера, а колеблются в окрестности исходного положения. Жидкость при этом в целом «макронеподвиж-
Рис. 1. Гармонические колебания трубы в цилиндрическом канале, заполненном вязкой несжимаемой жидкостью
на», т.е. среднее значение вектора скорости любой ее частицы стремиться к нулю при увеличении промежутка времени.
Для исключения времени используем комплексное представление перемещений и давления:
и = и о(г ,ф)е ; р = р о(г ,ф)е ,
(2)
где ш — круговая частота гармонических колебаний.
Подстановка выражений (2) в (1) дает
-рш2 и 0 =^р 0 + /шцУ2 и 0; V-и о = 0.
(3)
После перехода к проекциям в полярной системе координат уравнения (3) преобразуются к следующему виду:
-рш2ыг
-Рш иф =
дЫ„ и„
дР о дг
дР о гдф
диф
+ /шц + /шц
д2
д
д
2 \
дг 2 + гдг + г 2дф2
\
/д2
и;
^ г ?
дд + — +
2 \
дг2 гдг г 2дф2
иф; (4)
+ —+ — = о,
д г г гд ф
где иг, иф — проекции вектора перемещений и0 на радиальное и окружное направления полярной системы координат, т.е. амплитуды радиального и окружного перемещений.
ВШзехэшш] вьасзшш. ©аведжшй
Для исключения окружной координаты ф с учетом круговой геометрии и линейности задачи представим все искомые величины в виде одной гармоники разложения в тригонометрический ряд (первая гармоника): ur = ur 1(r )cos ф;
u ф = u ф1(г )sin ф; (5)
P 0 = Pi(r )cos ф.
Здесь индекс «1» соответствует первой гармонике. Обоснованием (5) служит то, что жесткое смещение внутреннего цилиндра в полярных координатах описывается аналогичным образом, т.е. в виде первой гармоники разложения в тригонометрический ряд: ur = A cos ф, uф = —A sinф. Из (4) и (5) следует система обыкновенных дифференциальных уравнений для компонент перемещений и давления:
-рю ur 1
dPi dr
+ /юц
1d 2ur 1 + dur 1
dr2
2 , P1 . •
-рю uф1 = +---h /юц
ф1 r
1d 2u
ф1
rdr du
\
\
dr2
+ ■
ф1
rdr
(6)
dur1 + + = 0
(7)
dr
Полученная система дополняется граничными условиями, которые имеют смысл условий прилипания жидкости к наружному и внутреннему цилиндру:
иг !</!) = А;
иг 1<г2) = 0;
и ф1(г1)=-а; и ф1(г2 ) = 0
Таким образом, для определения перемещений частиц жидкости и распределения давления в зазоре достаточно проинтегрировать систему (6) с граничными условиями (7), т.е. задача является полностью поставленной. Однако, система (6) не имеет аналитических решений. Численное решение этой системы (6) с помощью стандартных программ связано со следующими трудностями. Во-первых, старшие производные в уравнениях (6) умножаются на малый множитель ц, что обычно приводит к
численным проблемам. Во-вторых, структура (6) не позволяет непосредственно применять готовое программное обеспечение, что связано с кажущимся несовпадением порядка системы с количеством граничных условий. Если формально подсчитать сумму порядков старших производных в уравнение (6) с учетом dpl/dr, то получим тригонометрический ряд 5. Однако порядок системы равен 4 (в точном соответствии с количеством граничных условий), но программа решения систем дифференциальных уравнений обнаружить это не может и выдает сообщение об ошибке.
Таким образом, для численного решения систему (6) необходимо преобразовать. Для этого все неизвестные в системе (6) выражаются через одну неизвестную иг1(г) — амплитуду радиального перемещения: dur 1
ф1
■■—u„, — r-
dr
p1 = /юц
2 d ur, d ur, dur, r2-rL + 4r-r1 + -
\
\
dr3
dr2
dr
+
(8)
+рю2
2 dur 1 , \
r —— + rur 1 I.
\
dr
После исключения р1 и иф1 система (6) сводится к одному дифференциальному уравнению с одним неизвестным:
d4 иг1 _ /rd Зиг 1 ( 4 /гар'Ы2-
+ 6-3г +
rdr3
ц ) dr2
+
Ътр\ dur, 2ur,
1 r ■ + —f1 = 0.
(9)
\
цг ) dr
Порядок полученного разрешающего уравнения равен 4, следовательно, исходная система (6) имеет 4-й порядок.
Граничные условия также были представлены через радиальное перемещение:
иг 1 (Г1) = А;
иг 1 (Г2) = 0;
dur 1
dr
dur 1
dr
= 0;
= 0.
(10)
Миг1)
А
а
Рис. 2. Действительная (а) и мнимая (б)
Уравнение (9) с граничными условиями (10) адекватно воспринимается компьютерными математическими пакетами, так как порядок уравнения и количество граничных условий теперь одинаковы. Для решения линейной краевой задачи использована процедура NDSolve пакета МаШешаИса. Наличие малого параметра ц в знаменателях коэффициентов уравнения (9) делают уравнение «жестким», т.е. приводит к численным проблемам. Для решения проблемы жесткости применялся самый простой вычислительный прием, не требующий значительных изменений в алгоритме расчета — была значительно повышена точность всех
1т(иг1)
составляющие радиального перемещения
численных операций до 60 — 80 десятичных знаков. Пакет МаШешаИса позволяет выполнять вычисления и с большей точностью, но эффективность (быстрота) при этом значительно падает. Шаг интегрирования подбирается процедурой NDSolve автоматически, при этом вблизи границ интервала интегрирования (т.е. в зоне пограничных слоев) шаг интегрирования оказывается существенно меньшим, чем в других частях интервала интегрирования.
Типичные графики составляющих радиального и окружного перемещения представлены на рис. 2 и 3. Результаты получены при следующих значениях исходных данных, которые
Рис. 3. Действительная (а) и мнимая (б) составляющие окружного перемещения
ВШзехэшш] выкэдшпх ©аведжшй
примерно соответствуют параметрам тепловыделяющей сборки ядерного реактора ВВЭР-440 из работы [2]: г1 = 0,0745 м, г2 = 0,08 м, ц = =0,001 Пас, ш = 25 рад/с. Отличие заключается в геометрии канала. В [2] рассматривался канал шестиугольного сечения, а в данной работе канал имеет круглое сечение, поэтому соответствие только приближенное.
На рисунках 2, 3 хорошо видны пограничные слои. Кроме того, отчетливо видно, что наибольшее значение имеет действительная составляющая окружного перемещения, которую можно оценить из простейших геометрических соображений, использованных в работе [2]. Площадь сечения жидкости, выдавливаемая (ометаемая) одной четвертью внутреннего цилиндра радиусом г1 протекает через зазор Н = г2 — г1, поэтому амплитуда среднего по зазору относительного окружного перемещения должна быть равна <Ке(иф1)>/Л = = г1/Н = 0,08/0,0055 = 14,5. Примерно такое значение наблюдается в средней части графика Яе(мф1) на рис. 3. Заметна также небольшая асимметрия всех графиков относительно середины зазора, что и отличает рассматриваемое точное решение от приближенного из работы [2], где все решения либо строго симметричны, либо строго кососимметричны относительно середины зазора. При увеличении относительной ширины зазора отмеченная асимметрия значительно возрастает.
Коэффициент присоединенной массы и коэффициент демпфирования рассчитывались по формулам (11), аналогично работе [2]:
Яе[ Р1 (Г1
а = —-
$ =
рпг1 Аш 1т[ Р1 (г1 )]пг1 Аш
= 14,81;
(11)
= 362,2
Н-с
м
Приближенные формулы из работы [1] дают очень близкое к (11) значение для а и несколько отличающееся значение для $:
а = 14,88;
Н-с
$ = 339,4
м
Зависимости коэффициента присоединенной массы и коэффициента демпфирования от вязкости, найденные по предлагаемой методике, показаны на рис. 4, 5.
Крайние правые точки на рис. 4 и 5 соответствуют вязкости в 100 раз большей, чем вязкость воды. Большие вязкости, рассмотренные здесь, на практике встречаются очень редко. Однако следует учитывать, что цель применения результатов на практике для авторов статьи была не единственной. Авторы стремились также к получению максимально точных численных решений, на которых далее планируется тестировать конечно-элементные методики решения аналогичных задач. С этой точки зре-
Рис. 4. Зависимость коэффициента присоединенной массы от вязкости
Рис. 5. Зависимость коэффициента демпфирования
от вязкости
ния, полученные результаты вполне удовлетворительны, так как приближенные методики, основанные на рассмотрении пограничных слоев у каждого из цилиндров при больших вязкостях уже не применимы. Например, приближенная методика из работы [1] при вязкости жидкости ц = 0,1 Пас и тех же размерах и частоте дает значение а, завышенное почти в 2 раза и значение ^ заниженное в 3 раза, что при тестировании, очевидно, неприемлемо. В то же время при малых вязкостях приближенная методика из [1] дает практически точные значения а и что подтверждается и результатами авторов данной статьи.
Выводы
1. Разработана методика численно точного решения задачи о колебаниях длинной трубы в цилиндрическом канале, заполненном вязкой несжимаемой жидкостью.
2. Приведены результаты численных расчетов коэффициента присоединенной массы и коэффициента демпфирования для парамет-
ров, близких к параметрам тепловыделяющей сборки реактора ВВЭР-440.
3. Полученные результаты являются численно точными, поэтому они имеют не только практическое значение, но и являются хорошими тестами для отладки конечно-элементных методик.
Литература
1. Синявский В.Ф., В.С. Федотовский, А.Б. Кухтин, Л.В. Те-реник. Инерционность и гидродинамическое демпфирование при колебаниях труб и трубных пучков жидкости // Динамические характеристики и колебания элементов энергетического оборудования: Сб. ст. М.: Наука, 1980. С. 86—97.
2. Солонин В.И., Перевезенцев В.В., Сорокин Ф.Д. Демпфирование колебаний пучка твэлов чехловых тепловыделяющих сборок водоохлаждаемых реакторов// Вестник МГТУ им. Н.Э.Баумана. 2008. № 3. С. 75-85.
3. Шмелев В.Д., Драгунов Ю.Г., Денисов В.П., Василь-ченко И.Н. Активные зоны ВВЭР для атомных электростанций. Москва: ИКЦ «Академкнига», 2004, 220 с.
Статья поступила в редакцию 27.08.2012
