ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕУСТОЙЧИВЫХ РЕЖИМОВ РАБОТЫ ПРЕДОХРАНИТЕЛЬНОГО КЛАПАНА Текст научной статьи по специальности «Механика и машиностроение»

2020 Математика и механика № 68

УДК 533, 004

DOI 10.17223/19988621/68/13

Т. Редер, В.А. Тененев, А.А. Чернова

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

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

Ключевые слова: предохранительный клапан, численное моделирование, динамический процесс, метод С.К. Годунова, схема MUSCL.

Моделирование рабочих процессов в регуляторах давления на основе вычислительной динамики жидкости CFD (Computational Fluid Dynamics) дает основу для проектирования и оптимизации работы предохранительных клапанов. Основными затруднениями при осуществлении численных расчетов являются большое время вычислений и дискретное представление расчетной области при наличии подвижных твердых границ. В закрытом клапане имеются две изолированные области с условными названиями: нижняя - от входа до сечения седла клапана; верхняя - от седла до выхода. При подъеме диска с седла эти области соединяются через образовавшийся зазор. Возникает проблема выбора начальной величины зазора и соответствующей ему скорости движения диска. Уменьшение зазора приводит к соответствующему уменьшению расчетного шага по времени и к возрастанию времени счета, а при большом начальном зазоре теряется информация о предыстории движения диска. В разных работах проводятся разные обоснования начальной величины дискового зазора. Beune (2009) [1] проводил стационарные расчеты с зазором от 0.1 мм в 3D постановке. Динамика движения рассчитывалась в 2D-постановке, начиная с зазора в 0.01 мм. В работе [2] проводили SD-моделирование предохранительного клапана с применением ANSYS CFX. Начальный зазор составлял 0.1 мм. В этой работе подтверждено, что комбинированная модель способна предсказать динамику работы клапана, а также характеристики потока в разные моменты времени. Для расчета применялся метод динамического наложения слоев, требующий нескольких слоев сетки под диском. При численном моделировании работы клапана большинство авторов применяют пакеты ANSYS FLUENT и ANSYS CFX, а величина начального зазора меняется от 0.1 до 1 мм.

Поведение предохранительного клапана во время переходных режимов его работы (переходного процесса) является объектом исследования многих авто-

ров. В работе [3] анализируется математическая модель подпружиненного предохранительного клапана, соединенного с резервуаром сжимаемого газа через трубу. Клапан моделируется с использованием уравнений Ньютона и уравнений одномерной газодинамики. Одномерные модели с использованием эмпирических зависимостей гидравлической силы от характеристик течения применялись для исследования неустойчивых режимов работы клапана для жидкостных сред [4, 5].

Математическая модель работы клапана

Схема предохранительного клапана и расчетная область показаны на рис. 1.

Рис. 1. Расчетная область в клапане Fig. 1. Valve computational domain

Расчетная область разделена на две части. Первая часть является осесиммет-ричной (l = 1) и в ней используется цилиндрическая система координат (х, r, ф). Во второй части (l = 0) применяется прямоугольная система координат (х, y, z). Безразмерные уравнения, описывающие движение вязкого теплопроводного газа, имеют следующий вид:

d(ql2m) s[q2 (A - Av )] э[q2 (B - Bv )] 5(c - Cv )

- = S.

(1)

дt дц1 дд2 дд3

В зависимости от принадлежности к части расчетной области, координатами Ч = (?1,Ч2,Ч3) являются:

[(х, у, z) 11 = 0

К г, ф)|/ = 1,

где величина I = 0 соответствует декартовой системе координат (при этом У = 1), а I = 1 - цилиндрической (г1 = г);

(1, q2, qs ) = <[(

т

р\ ри 1 { ру 1 рw

ри рии + р риу рuw

ру руи руу + р руw

р^ , А = рwu , Б = рwу , С = рww + р

е и(е + р) у(е + р) w(e + р)

рк риК руК рwK

ре V риЕ ч руЕ у рwE

А =-

Яе0

121

Р. у + л±*1рТ|1

стт дq1

Г X 1 дК

А+-

дql дЕ

дq1

Б,

Яе0

12

п22

Р2. V+ л *

стт ¿^ ^ X 1 дК

Л+-

дq2 дЕ

дq2

С„ =

Яе0

Чэ

Т23

1зз

стт дq3

л+-

5qз дЕ 5qз

8 =

О О

/1 р + рт^2 -

33

Яе0

11 -ру^+ —23 ' Яео

q2р(2 - Е)

q2 Кр(-с2е)

Вектор скорости V = (и,у, ^)т имеет составляющие в соответствующем базисе (q1,q2,q3). Векторы т, А,Б,С, Ау,Бу,Су,8,8у содержат газодинамические

р у2

комплексы, составленные из переменных: р - плотность газа; е = ре + р1

2

полная энергия единицы объема газа; внутренняя энергия определяется уравнени-

гт Р СР ем состояния совершенного газа е = су1 = --—, где у =--показатель адиа-

(у-1)р Су

баты; срсу - изобарные и изохорные теплоемкости; Т - температура, А, Б, С, -конвективные члены: Ау,Бу,Су - вязкие члены; 8 - источниковый член.

Переменные р, V отнесены к критическим значениям р*, а*, давление p и

2 I Р*

энергия e к р*а* , где а* = /у— - скорость звука. Пространственный масштаб

р*

r

r0 - радиус седла клапана, масштаб времени т = —. Критические значения пере-

a*

менных связаны с параметрами торможения в контролируемой емкости известными соотношениями:

Y 1

( Y + 1 | Y-1 ( Y + 1 | Y-1

Р* =Ро J , Р* = Ро|

В качестве масштабов энергии турбулентности K и скорости ее диссипации E

выбраны — и —2 соответственно. Турбулентная вязкость, отнесенная к v0 : т т2

X = c K—. Коэффициент кинематической вязкости —0 соответствует температуре

E

торможения T0. Зависимость вязкости от температуры учитывается формулой Саттерленда в виде

L/ о-5 Y + 1 | y- 1 ( e | e0 + e,

— (e) = —0Л, Л = —0 (l+J^)

e + e„

где в0 = ^ в, = ^, Т = 122 К.

а* а*

Компоненты П] образуют тензор

П = р(/ + 2ёе/ (V)-|<Ну (V)), Б = йе/ (V),

Рг =К ) , Пе = С = С1 1^1 -Ле.

1 Е 11 1 + вЕЦ1 { Пе о

Константы имеют следующие значения:

сц-0.085; с1 = 1.42; с2 = 1.68; РЕ = 0.012; Пео = 4.38; ак = 0.72; а, = 0.72,

т-т ^ г0 а* „ „ аТ - число Прандтля; Яе0 =--число Реинольдса.

Перемещение диска в осевом направлении , совпадающим с осью х определяется действием силы со стороны газа Г/; силы упругости пружины Г, = К, ); демпфирующей силы трения = Кёп и силы тяжести .

Здесь: К, - коэффициент жесткости пружины; - перемещение диска (высота подъема); - начальное сжатие пружины: Кё - коэффициент демпфирования (трения); п - скорость перемещения диска.

Уравнения движения диска имеют вид

ё п

ш„

dt

= Ff - Fs- Fd- m

dt

= П

(2)

(3)

Адиабатическое изменение давления торможения в емкости, откуда происходи,

дит сброс газа, определяется изменением массы газа в емкости

dt

- = Gin - G

valve

при расходе газа через клапан Ота1уе и приходе газа от внешних источников О1Г

dp =]RTnr -G )

dt W K In valveJ l p0

Y± P I Y

(4)

где Ж - объем емкости; Я - газовая постоянная; Т0, р0 - температура и давление газа в начальный момент времени.

Начальные условия для решения уравнений газодинамики соответствуют закрытому клапану. Это значит, что в объеме под диском клапана параметры равны параметрам в емкости, а в объеме над диском - параметрам окружающей среды.

Граничные условия для уравнения (1):

- Вход. Задаются направление вектора скорости, условие постоянства энтропии, соотношения на характеристике, входящей в область. Для К и Е заданы начальные уровни турбулентности, соответствующие х = 1.

- Выход. Задано давление ра, для параметров турбулентности ставятся «мягкие» условия.

- Твердые стенки. Для скорости - условие прилипания, для К и Е - нулевые градиенты по направлению нормали.

Начальные условия: до зазора давление, как в резервуаре; после зазора - давление окружающей среды; скорости нулевые.

Для численного решения систем газодинамических уравнений (1), применяется метод контрольного объема. Для каждой грани с индексом И выражение для потоков имеет вид

" Я (у-Ы) Я (у-Ы )и-пхР Я (у-Ы)У-пуР

Я(у-Ы)Ж-п2Р , (5)

Я (у-Ы )О-ЫР Я (у-Ы )У Я (у-Ы)

где Я,и,У,Ж,О,Р,У,2 - величины, соответствующие плотности, скоростям, энергии, давлению, энергии турбулентности, скорости диссипации; у - скорость грани в направлении нормали; Ы - нормальная составляющая вектора скорости; стИ - площадь поверхности соответствующей грани.

Qh =At4

Для каждой грани решается соответствующая задача о распаде разрыва с двумя наборами газодинамических параметров в контрольных объемах, разделенных этой гранью. Параметры газа на границах контрольных объемов определяются по методу С.К. Годунова с использованием автомодельного решения задачи о распаде произвольного разрыва [6]. К сожалению, способы приближенного решения задачи Римана, сокращающие время расчета [7], хорошо работают только в относительно простых областях и для тестовых задач. В случае осесимметричного течения в клапане появляются осцилляции в давлении. В трехмерной постановке для условий сложного течения в клапане схема становится неустойчивой. Для повышения порядка аппроксимации разностного метода Годунова применяется схема Ми8СЬ. Значения газодинамических параметров для решения задачи о распаде разрыва определяются с применением экстраполяции с ограничителем [8, 9]. Например, для направления по индексу г для правой грани:

- значения переменных / = (р, и, V, м>, р) «справа» от поверхности разрыва

/к = /г +1 к + 0.5Т(С)(/1]к - /+1 ]к ) ,

— /г+2 ]к

С =

fijk fi+1 jk

Y(c ) =

2

С + С

1+С2

- значения переменных «слева» от поверхности разрыва

fL = j + 0.5T(c )(f+i ]к - j ) , c = ffi .

f+1 jk Jijk

Величины R,U,V, W,G определяются из решения общей задачи о распаде произвольного разрыва. Величины Y, Z на гранях контрольного объема определяются с применением экстраполяции с ограничителем в соответствии с направлением нормали скорости.

Для уравнений (1), записанных в цилиндрической системе координат, строится разностная сетка в блоке 1 в плоскости ф = const комплексным методом граничных элементов. Контрольные объемы представляют собой секторальные вырезки из цилиндров с криволинейной образующей (рис. 1). Для таких контрольных объемов нормаль и касательная к боковым граням рассчитываются просто, так как сетка является ортогональной. Для блока 2 контрольные объемы могут иметь вид произвольного многогранника, причем грани не обязательно плоскости. Для пространственного четырехугольника с вершинами 1, 2, 3, 4, пронумерованными против часовой стрелки, вектор нормали имеет составляющие:

"AV21Az32 - AV32Az21 + A>'43Az14 - ^14^43 "

Az21A32 -Az32A21 + Az43A14 -A14A43

_ЛГ21ДУ32 -Ax32Ay21 + Ax43ДУ14

где, например, Ax21 = x2 - x1. Тангенциальный вектор возьмем в виде

nx

x

n = ny = 0.5

_ nz _

T x Ax31

T = T y = AJ31

_Az31 _

Векторы п, Т нормируются и третий орт определим как векторное произведе-

ние

К =

= [ х Т].

Проекции вектора скорости на орты, необходимые для решения задачи о распаде разрыва, это скалярные произведения N = (V, п); Т = (У,Т); К = (У,К).

После решения задачи о распаде произвольного разрыва восстанавливаются величины скорости на ребре:

и = п^ + тХТ + кхК , V = п^ + туТ + куК , Ж = пN +хгТ + к2К

и вычисляются потоки (5) через соответствующие грани.

Градиенты переменных, входящие в компоненты , тензор П вычисляются

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

Л-=-а ъ,

где

Ъ1 = Л+1 ]к Л]к , С1 = Чг+1 ]к %'к , Ъ2 = /¡]+1к -1к + Л+1 ]+1к Л+]-1к , С2 = Чг]+1к - Чг]-1к + Чг+1 ]+1к - Чг+1 ]-1к , Ъ3 = Л]к+1 - Л]к-1 + /г+1 ]к+1 - /г+1 ]к-1 '

С3 = Чг]к+1 Чг]к-1 + Чг+1 ]к+1 Чг+1 ]к-1 , С = С1 ' (С2 Х С3 ) , С =

у

С

Переход к следующему шагу по времени осуществляется по двухшаговой схеме со вторым порядком точности:

1 1

п+— п+—

щ.] = т п

Г^п I \ ' г\п _ о п гп

г,],кОг,],к -1 I Ой - 8г,],кОг,],к

АГп

п+1 п+1 _ тг, ] ,кОг, ] ,к = т

п

г,],кОг,],к '

1 1 П

__ п+— п+— п+—

I Ог 2 -

2

АГп

(блок 1)

1

п+—

тг, Д = т"] ,к"

О

1 ÍS ог - К]*}^,

т

п+1

г, ],к

т

1

г,],к V й (

г,],к '

о,

г, ],к

,] 2 1 Л

I Ог 2 - 8 2

п+-г, ],к

АГп

/

(блок 2)

где, г,],к - номер контрольного объема; О" ^ к - величина контрольного объема на временном шаге п ; суммирование потоков проводится по всем граням контрольного объема. В первом блоке контрольные объемы деформируются при движении диска, во втором в течение времени не изменяются.

Задача совместного решения уравнений (1) - (3) является сопряженной. Связь осуществляется через газодинамическую силу, определяемую через интеграл давления и вязкого напряжения трения по поверхности диска

р/ = {(Уйгй ф

и деформацию подвижной границы блока 1. Величина вязкого напряжения трения т№ определяется через пристеночные функции [11].

Дифференциальные уравнения (2),(3) интегрируются по неявной схеме [2]

( £ п — £ П-1 £ п \

(( - К £о - т^ ) + К,е £п + т:, + А I

£ т

к А/ + -А- + К, А/

Анализ начального участка движения диска

Цикл работы предохранительного пружинного клапана начинается от момента времени, когда давление в контролируемом резервуаре превышает допустимое. Газодинамическая сила, действующая на запорный диск клапана, по величине превышает силу упругости пружины. Диск начинает перемещаться, площадь проходного сечения увеличивается и рабочая среда из резервуара сбрасывается через выходной патрубок в атмосферу. Циклограммы для клапана 213 «давление - перемещение диска» показаны на рис. 2.

10.8

■---experiment

■ simulation

10.4

10

4

9.6 -

9.2

4 6

Перемещение диска, мм

2

3

1

Рис. 2. Циклограмма клапана «давление - перемещение диска» Fig. 2. «Pressure - disc displacement» cyclogram of the valve

Из точки 1 давление возрастает до предельного значения в точке 2 и клапан начинает открываться. Точка 3 соответствует максимальному ходу диска. На этом рисунке принято избыточное давление. В точке 2 преодолевается сила упругости пружины. С увеличением зазора между диском и седлом клапана газодинамическая сила возрастает, все время превышая силу упругости пружины. В точке 3 диск занимает предельное положение и начинается более интенсивный сброс давления в резервуаре, приводящий к снижению газодинамической силы. При достижении газодинамической силы величины силы упругости пружины в точке 4 диск начинает опускаться. Газодинамическая сила снижается как за счет снижения давления, так и за счет уменьшения дискового зазора. Цикл замыкается с закрытием клапана и приходом в точку 1.

Газодинамические процессы, требующие численного моделирования, соответствуют участку циклограммы 2-3-4. При проведении численного моделирования динамических процессов в пружинных клапанах и анализе результатов все авторы отмечают проблематичность расчета стартового отрезка времени. При величине зазора порядка 0.01 мм при количестве десяти слоев дискретный шаг по времени порядка 10-8 5, в зависимости от допустимого числа Куранта, что приводит к большим временным затратам даже на высокопроизводительных компьютерах.

Запишем уравнение (2) в виде:

= ^Р) - К (^ + §0) - ксП - . (6)

Начальные условия для уравнений (2), (3): Г (0) = 0 ; п(0) = 0.

Предохранительный клапан настраивается на срабатывание при достижении заданного уровня давления Р0 . Сила со стороны газа на диск в момент достижения давления Р0 равна

^ (Г, Р) = ст(Р0 - Ра), (7)

где ра - противодавление окружающей среды; ст - площадь проходного сечения

седла клапана, определяемая конструкцией клапана. При превышении газодинамической силы упругости пружины и силы веса подвижной части ст(р0 - ра) > К550 + т^ начинается движение диска в соответствии с уравнением (6).

Разложим газодинамическую силу в ряд Тейлора в окрестности точки (* = 0;0;р = Р0):

Р (Г Р)= Р (0 Р ) , дЕЛ (ГР) Г . дЕГ (ГР) СР| * + 0 (*2 Ч

(г,Р ) = (0, Р0 )+-дг-1г=0, Р= Р0 г+-дР-г=0, Р= Р0 ~^\*=0* + 0 (г , * ).

Из выражения (7) следует

дР; (Г, Р)

дР

Г=0, Р=Р0 = ст,

ф| ^ кЯТ0

а из (4) -СЦ = От-Ж1

т Ж

с учетом того, что Ота1те (0) = 0 .

Производную газодинамической силы по аргументу заменим конечно-разностным соотношением:

^ (£, р)■ ^ ¥/ (А£, Ро)- (0, Ро) = А¥

^=0, р=Ро ~ А§ А§ ,

где - некоторая величина зазора.

Результаты анализа численных расчетов и экспериментов [12] показывают, что за время подъема диска на максимальную величину давление в емкости практически не изменяется. Поэтому при анализе срабатывания клапана важной является

величина производной —к=0 , обеспечивающая выход динамической системы,

Ж

описываемой уравнениями (2) - (4), из состояния покоя. Величина силы (А|, р0) при заданной величине зазора вычисляется при численном решении газодинамической задачи.

С учетом того, что (0, р0) = К,50 + ш^ , А¥ = (А§, р0) - (0, р0), для

перемещения диска получаем уравнение

й\ в ,, —+ к—= а£+Ь/, С2 сН

где введены обозначения:

к = ^, а = ГА¥ - к 1-1, Ь = ШТ0°п

у ш,

Решение линейного неоднородного обыкновенного дифференциального уравнения второго порядка имеет вид

£ = Се^ + С^2 + , где Х12 = -±Л/— + а .

1 2 а 1,2 2 V 4

Из начальных условий (0) = 0 , ) = 0 определяем константы интегриро-

&

вания С1, С2 и получаем аналитические выражения для зависимостей перемещения и скорости диска от времени:

" -К ил)

а

(8)

4а

П = + (9)

, к2 где А =--+ а .

4

Из уравнения (8) определяется время достижения заданной величины зазора , из уравнения (9) - скорость движения диска в этот момент времени. Для дальнейшего решения эти величины применяются в качестве начальных условий при совместном решении уравнений (1) - (3).

Результаты расчетов полного цикла работы клапана

При движении диска сетка перестраивается. Для деформации сетки в процессе численного расчета предварительно рассчитываются несколько вариантов сетки от начального до конечного значения зазора (для расчета взято 10 положений диска). При моделировании по текущей величине зазора строится текущая сетка по двум промежуточным вариантам с помощью интерполирования. Сетка во втором блоке является неструктурированной и не деформируется. Построенная таким образом сетка отслеживает участки течения с большими градиентами. Кроме того, такая сетка, близкая к ортогональной, легко измельчается или укрупняется при необходимости.

С применением рассмотренного алгоритма проведены расчеты полного цикла работы клапана 213 при максимальном абсолютном давлении в резервуаре 11.6 бар. На рис. 2 показаны расчетные и экспериментальные циклограммы «давление -перемещение диска». Условия проведения экспериментов описаны в работе авторов [12]. Сравнение расчетных и экспериментальных циклограмм «давление - перемещение диска» является комплексным показателем возможности моделирования полного цикла работы клапана. В рассчитываемых условиях) клапан работает в режиме полного подъема диска в течение основного времени процесса. Характеристикой конструкции клапана для различных давлений служит зависимость приведенной силы (отношение газодинамической силы к величине давления) от высоты подъема диска. Результаты численного моделирования и экспериментальных измерений представлены на рис. 3.

180

160 -

я а с й ч я О

140 -

120 -

100

-simulation

О experiment

2 4 6 8

Подъем, мм

Рис. 3. Зависимость приведенной силы от высоты подъема диска Fig. 3. Dependence of the reduced force on the height of the disc

10

Из рис. 3 следует, что расчетные значения приведенной силы на основном участке работы находятся в пределах доверительного интервала при доверительной вероятности 0.99. При малой величине подъема диска наблюдаются большие разбросы экспериментальных значений. Эти разбросы связаны с вибрациями при по-

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

Расчеты показывают сложную структуру течения (рис. 4).

Рис. 4. Структура течения и турбулентная вязкость в клапане Fig. 4. Flow structure and turbulent viscosity in the valve

В левой части показаны линии тока и величины числа Маха. Течение после зазора становится сверхзвуковым с образованием ударных волн. Большая деформация поля скорости приводит к росту турбулентной вязкости (правая часть рис. 4).

Нестабильные режимы работы клапана

При наличии подводящей трубы от резервуара к клапану в экспериментах наблюдается колебательный режим движения диска и давления в клапане. Одномерный расчет таких режимов проводился рядом авторов [3-5] при различных допущениях о связях изменения газодинамики с движением диска. Трехмерный газодинамический расчет на сетке порядка 106 узлов (рис. 1) с шагом 10-8 с на полный цикл работы 3 с является очень трудоемким даже при распараллеливании вычислений. Для возможности оперативного проведения расчетов динамического режима работы клапана предлагается перейти от трехмерной вязкой постановки к двумерной невязкой в блоке 1 на грубой сетке. Коррекция полученного грубого решения осуществляется введением корректирующих множителей для газодинамической силы и массового расхода, полученных из сопоставления с трехмерным расчетом. Грубая сетка строится из исходной в блоке 1 отбрасыванием коорди-

натных линий. Так как продольные координатные линии полностью отслеживают контур диска, то возможо рассчитать газодинамическую силу, хотя течение здесь следует считать промежуточным между одномерным и двумерным. Из сопоставления с трехмерным вариантом расчета определяются:

1 - поправочная функция к нижней приведенной силе, зависящая от перемещения диска |: у(§) = 1 + 0.0156^+ 0.0009|2;

2 - поправочная функция к массовому расходу g (|) = 1 + 0.0584^ - 0.0045|2;

3 - поправочная функция к верхней силе / = 1 - 0.00811 .

Полученные в расчете 2Б значения параметров умножаются на соответствующие поправочные функции. При данном алгоритме трехмерность течения учитывается поправочными функциями, а нестационарность обеспечивается совместным решений уравнений (2), (3) и нестационарных 2Б-уравнений газодинамики на грубой сетке.

Для повышения производительности также применялись параллельные вычисления с разделением расчетной области 1 с осевой симметрией на количество блоков в продольном направлении, равное количеству потоков (логических процессоров). Границы блоков являются динамическими для обеспечения равенства времени расчета во всех блоках.

Влияние длины подводящей трубы на стабильность работы клапана

Рассмотрим влияние длины подводящей трубы на колебательные процессы в клапане. Длина трубы изменялась от 0 до 2 м. Исследуется режим, когда диск не прижимается полностью к верхнему ограничителю. Движение диска показано на рис. 5.

12

tfl

«

о

С

Рис. 5. Зависимость перемещения диска от времени Fig. 5. Disc displacement as a function of time

Номера кривых соответствует длине трубы: 1 - Ь = 0; 2 - Ь = 0.5 м; 3 -Ь = 0.75 м; 4 - Ь = 1 м; 5 - Ь = 1.5 м; 6 - Ь = 1.75 м. С увеличением длины трубы амплитуда колебаний диска возрастает, а частота колебаний уменьшается (рис. 6).

8

4

А, мм

0

Т, мс 18

16

14

12

0

0.5

1.5

L, м

0.5

1

1.5

L, м

Рис. 6. Амплитуда колебаний диска (а); величина периода колебаний (b) Fig. 6. (a) Disc oscillation amplitude and (b) the oscillation period

a

2

1

1

b

Причиной колебаний является образование волн давления в трубе после открытия клапана (рис. 7). Показано изменение давления на входе в клапан и величины газодинамической силы.

Собственная круговая частота колебаний подвижной части клапана равна 154.7 1/с с длиной периода Т = 40.6 мс большей в два раза, чем на рис. 6, Ь. Для определения характера движения газа в трубе проведен расчет при фиксированном зазоре 1 мм. Начальные условия такие же, как при открывании клапана. Изменение давления во времени в различных сечениях трубы показано на рис. 7, Ь. Волна доходит до входа в трубу, где со временем давление стабилизируется. На входе в клапан давление также стабилизируется. В середине трубы происходит волновой процесс, связанный с отражением волн. Период колебаний примерно 10 мс, что примерно в два раза меньше, чем при движении диска. Частота колебаний при движении диска определяется сложным характером взаимодействия звуковых волн с упругой системой диска. В короткой трубе колебания быстро затухают и процесс стабилизируется. При длине трубы 2 м клапан начинает открываться, а затем на фазе снижения давления в результате волнового процесса полностью закрывается. Далее он может открыться, если давление в резервуаре возрастет.

40 30

20 10

0 20 40 60 80 t, мс

p, бар 22.2

21.8 21.4 21.0

20.6

0 5 10 15 20 25 t, мс

Рис. 7. Волновой процесс в клапане для длины трубы 1.75 м, зависимости от времени: а - перемещения диска, мм (кр. 1); давления, бар (кр. 2); газодинамической силы, кН (кр. 3); силы пружины, кН (кр. 4); b - давления, бар: на входе в трубу (кр. 1); в середине трубы (кр. 2); в конце трубы (вход в клапан) (кр. 3) Fig. 7. Wave process in the valve for a pipe of 1.75 m length. (a): 1, disc displacement in mm as a function of time; 2, pressure in bar; 3, gas dynamic force in KN; and 4, spring force in KN; (b): 1, pipe inlet; 2, middle of the pipe; and 3, pipe outlet (valve inlet)

Заключение

Для численного моделирования работы предохранительных клапанов применен метод С.К. Годунова в трехмерной постановке для вязкого теплопроводного газа. Для повышения точности с сохранением монотонности решения используется схема MUSCL (Monotone Upwind Schemes for Conservation Laws) и двухшаго-вый по времени метод. Сравнение численных расчетов газодинамической силы с результатами проведенных экспериментов для клапана 2J3 показало их удовлетворительное согласование. Существующая проблема расчета стартового отрезка времени при малой величине дискового зазора решается с применением аналитического решения, используемого для определения начальных условий при моделировании. По результатам численных 3Б-экспериментов динамических режимов работы клапана осуществлен переход к нестационарной 2Б-модели на грубой осесимметричной сетке, дающий согласование расчетной и экспериментальной

b 1

\ 2 1 \ / J3 Г

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

ЛИТЕРАТУРА

1. Beune A. Analysis of high-pressure safety valves. Eindhoven: Technische Universiteit Eindhoven, 2009. 134 p. DOI: 10.6100/IR652510.

2. SongX.G., Wang L.T., Park Y.C., Sun W.A. fluid-structure interaction analysis of the spring-loaded pressure safety valve during popping off // 14th International Conference on Pressure Vessel Technology. Procedia Engineering. 2015. V. 130. pp.87-94. DOI: 10.1016/j.proeng. 2015.12.178

3. Hos C.J., Champneys A.R., Paulc K., McNeelyc M. Dynamic behavior of direct spring loaded pressure relief valves in gas service: Model development, measurements and instability mechanisms // Journal of Loss Prevention in the Process Industries. 2014. V. 31. Р. 70-81.

4. Hos C.J., Champneys A.R., Paul K., McNeely M. Dynamic behaviour of direct spring loaded pressure relief valves: III valves in liquid service // Journal of Loss Prevention in the Process Industries. 2016. V. 43. P. 1-9. DOI: 10.1016/j.jlp.2016.03.030.

5. Dimitrov S., Komitovski M. Static and dynamic characteristics of direct operated pressure relief valves // Machine Design. 2013. V. 5. No. 2. P. 83-86.

6. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976. 400 c.

7. Сафронов А.В., Фомин Ю.В. Метод численного решения уравнений газодинамики с помощью соотношения на разрывах // Труды МФТИ. 2010. Т. 2. № 2. С. 137-148.

8. van Albada G.D., van Leer B., Roberts W.W.Jr. A comparative study of computational methods in cosmic gas dynamics // Astronomy and Astrophysics. 1982. V. 108. No. 1. P. 76-84.

9. Wesseling P. Principles of Computational Fluid Dynamics. Berlin: Springer-Verlag, 2001. 644 p.

10. Wesseling P., Segal A., and Kassel C.G.M. Computing flows on general three-dimensional nonsmooth staggered grids // Journal of Computational Physics. 1999. V. 149. P. 333-362.

11. Волков К.Н., Емельянов В.Н. Моделирование крупных вихрей в расчетах турбулентных течений. М.: Физматлит, 2008. 368 с.

12. Редер Т., Тененев В.А., Паклина Н.В. Исследование влияния величины начального зазора на динамику открывания предохранительного клапана // Интеллектуальные системы в производстве. 2018. Т.16. № 2. С. 28-40. DOI 10.22213/2410-9304-2018-2-28-40.

Статья поступила 07.10.2019

Raeder T., Tenenev V.A., Chernova A.A. (2020) NUMERICAL SIMULATION OF UNSTABLE MODES FOR A SAFETY VALVE. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 68. pp. 141-157

DOI 10.17223/19988621/68/13

Keywords: safety valve, numerical simulation, dynamic process, Godunov method, MUSCL scheme.

When designing pressure regulators, one needs to have a complete understanding of gas-dynamic processes. The numerical algorithm for three-dimensional gas-dynamic modeling of a full cycle of spring safety valve operation is proposed, which allows one to significantly reduce the computing time. Grid reconfiguration during CFD modeling is provided by interpolation

procedure using previously calculated grids. Calculations show that gas-dynamic numerical simulation should account for a three-dimensional structure of the unsteady flow and the motion of the disc. These factors are taken into account when calculating full cycle of the valve on a coarse grid with the use of correction functions for the force and flow characteristics of the valve. The correction functions are calculated by the false transient method in the three-dimensional formulation. Cyclograms of the valve operation demonstrate satisfactory agreement of the experimental and numerical simulation results. The agreement in the variation of gas-dynamic forces with time is observed, except for the transitional regime before the valve starts to close. In the main work area, the calculated values of the reduced force belong to a confidence interval.

Tomas RAEDER (Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation). E-mail: raeder.t@leser.com

Valentin A. TENENEV (Doctor of Physics and Mathematics, Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation). E-mail: v.tenenev@gmail.com

Alena A. CHERNOVA (Candidate of Technical Sciences, Kalashnikov Izhevsk State Technical University, Izhevsk, Russian Federation). E-mail: alicaaa@gmail.com

REFERENCES

1. Beune A. (2009) Analysis of High-pressure Safety Valves. Eindhoven: Technische Universiteit Eindhoven.

2. Song X.G., Wang L.T., Park Y.C., Sun W.A (2015) Fluid-structure interaction analysis of the spring-loaded pressure safety valve during popping off. Procedia Engineering. 130. pp. 8794. DOI: 10.1016/j.proeng.2015.12.178.

3. Hos C.J., Champneys A.R., Paulc K., McNeelyc M. (2014) Dynamic behavior of direct spring loaded pressure relief valves in gas service: Model development, measurements and instability mechanisms. Journal of Loss Prevention in the Process Industries. 31. pp. 70-81. DOI: 10.1016/j.jlp.2014.06.005.

4. Hos C.J., Champneys A.R., Paul K., McNeely M. (2016) Dynamic behaviour of direct spring loaded pressure relief valves: III valves in liquid service. Journal of Loss Prevention in the Process Industries. 43. pp. 1-9. DOI: 10.1016/j.jlp.2016.03.030.

5. Dimitrov S., Komitovski M. (2013) Static and dynamic characteristics of direct operated pressure relief valves. Machine Design. 5(2). pp. 83-86.

6. Godunov S.K., Zabrodin A.V., Ivanov M.Ya., Krayko A.N., Prokopov G.P. (1976) Chislennoe reshenie mnogomernykh zadach gazovoy dinamiki [Numerical solution of multidimensional gas dynamics problems]. Moscow: Nauka.

7. Safronov A.V., Fomin Yu.V. (2010) Metod chislennogo resheniya uravneniy gazodinamiki s pomoshch'yu sootnosheniy na razryvakh [Method for the numerical solution of gas dynamics equations using the discontinuity relation]. TrudyMFTIProceedings ofMIPT. 2(2). pp. 137-148.

8. van Albada G.D., van Leer B., Roberts W.W.Jr. (1982) A comparative study of computational methods in cosmic gas dynamics. Astronomy and Astrophysics. 108(1). pp. 76-84.

9. Wesseling P. (2001) Principles of Computational Fluid Dynamics. Berlin: Springer-Verlag.

10. Wesseling P., Segal A., Kassel C.G.M. (1999) Computing flows on general three-dimensional nonsmooth staggered grids. Journal of Computational Physics. 149. pp. 333-362. DOI: 10.1007/978-94-017-1564-5_2.

11. Volkov K.N., Emel'yanov V.N. (2008) Modelirovanie krupnykh vikhrey v raschetakh turbulentnykh techeniy [Modeling of large vortices in calculations of turbulent flows]. Moscow: FIZMATLIT.

12. Reader T., Tenenev V.A., Paklina N.V. (2018) Trekhmernoe chislennoe modelirovanie gazodinamiki predokhranitel'nogo klapana [Numerical 3D simulation of safety valve gas dynamics]. Vestnik IzhGTU imeni M.T. Kalashnikova Bulletin of Kalashnikov ISTU. 21(4). pp. 174-181. DOI: 10.22213/2413-1172-2018-4-174-181.

Received: October 7, 2019