CC BY
27
УДК 66.012-52
С. А. Мерзляков, Д. В. Елизаров, В. И. Елизаров, С. Г. Дьяконов
МОДЕЛИРОВАНИЕ СТРУКТУРЫ ПОТОКА ЖИДКОСТИ НА КОНТАКТНЫХ УСТРОЙСТВАХ СРЕДСТВАМИ ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИКИ
Ключевые слова: контактное устройство, уравнение сохранение импульса, вычислительная гидродинамика, 3D моделирование.
На основе программных систем геометрического моделирования, конструирования, средств вычислительной гидродинамики, математического моделирования рассматривается алгоритм определения структуры потока и моделирования ситчатого контактного устройства.
Keywords: contact device, the equation of conservation of momentum, computational fluid dynamics, 3D modeling.
Based on software systems of geometric modeling, engineering, computational fluid dynamics tools, mathematical modeling is considered an algorithm for determining the flow and structure modeling perforated contact device.
Разработка конструкции контактного устройства в настоящее время проводится методами физического, математического и гидродинамического моделирования. Наиболее перспективным методом разработки конструкции контактного устройства является интерактивный режим диалога конструктора и программных вычислительных комплексов, который можно поделить на два этапа.
Моделирование макета ситчатой тарелки в програмной среде Ansys
Для разработки модели контактного устройства проектировщик может использовать различные программные системы 3D моделирования, как внутренние компоненты среды Ansys - Design Modeler, так и внешние, например, Компас, AutoCad. Для гидродинамического моделирования структуры потоков жидкости и газа в барботажном слое контактного устройства используется среда Fluent системы Ansys.
Первоначально проектировщик составляет трехмерную твердотельную полноразмерную (3D) модель в среде моделирования Компас или AutoCad, Ansys. Для упрощения задачи моделирования контактного устройства в силу симметрии движения потоков для примера рассматривается модель половины ситчатой тарелки диаметром 700мм, ограниченная областью барботажного слоя без переливных устройств [1, 2].
Затем проект контактного устройства в виде файла формата (sat) передается в программный комплекс Ansys. В области барботажного слоя методом конечных элементов проводится построение сетки (Ansys Mesh), назначаются входные, выходные поверхности, стенки, начальные условия (количество ядер, точность расчета), выбирается тип решаемой задачи (стационарная, нестационарная), для барботажного слоя Multiphase Model (модель Эйлера), стандартная k — s модель, задаются граничные условия на стенке, оси симметрии, входных и выходных областях газа и жидкости, вводятся начальные значения для инициализации решения (скорость диссипации энергии, масштаб турбулентных пульсаций, кинетическая энергия). Расчет гидродинамических параметров ведется по нестационарным
уравнениям переноса: уравнение сохранения импульса в инерциальной системе отсчета:
^ (рб ) + У- (рб б) = -Ур + У-(Т| + рд + Р,
где р - статическое давление, т - тензор напряжений, рд и Р - гравитационные силы и внешние массовые силы, рб - импульс.
Тензор напряжений
т = |
(vu + VuT )-
(Vu + VuT )— 2V-uI
где | - молекулярная вязкость, I - единичный тензор. Стандартная k — s модель:
k )+^(Р kuj) =
at
д x
a x
i +
+ Gk + Gb — Ps— YM + Sl
dt
(ps)+aj (psuj)=
dx.
lt
. CTsy
I1 k
ds
dk
д x i
дХ,
pqss—pC2s
--/= + C1s7 C3sGb + Ss'
k+Vvs k
C1 = max
0.43,-
Л + 5
sk
л = S-,
s
где О1в = 1,44; О2в = 1,92; Оц = 0,09; стк = 1; стЕ = 1,3
- число Прандтля кинетической энергии турбулентности и диссипации, соответственно; Gk - производство кинетической энергии турбулентности, G ь
- производство кинетической энергии турбулентности за счет плавучести, О 3 в - степень влияния в на плавучесть.
Кинематическая турбулентная вязкость опреде-
ляется как
C k2
it =PC| —
Уравнение переноса для площади межфазной поверхности:
д(Р gXp ) „ / - \ 1 DPg
+v-(p gu g X p )= X p +
p g (S RC + S WE + S TI )
2m
3 ag
-X p +
д
д
д
+
2
Л
где Хр - площадь межфазной поверхности (т2/т3), ад - объемная концентрация газа, тд скорость массопереноса газовой фазы в единице
Таблица 2 - Начальные условия
объема смеси (kg/m3/s), SF
^с и SWE - члены уравнения отвечающие за теплоотдачу смешением за счет случайного столкновения и вследствие уноса, ST| источник распада.
Взаимодействие фаз описывается моделью Шил-лера-Наумана
f = CDRe
24
2411 + 0.15 Re0 687 I/Re , Re < 1000
l(i + 0.
I 0.44 Re > 1000
Критерий Рейнольдса Re для сплошной фазы q и дисперсной p:
Pq I" p
Re =-
"qdp
, q
При решении нестационарной задачи вводится время дискретизации, количество итераций на каждом шаге и количество шагов дискретизации. Процесс расчета считается законченным, если число итераций на последнем шаге дискретизации не превышает 5-10.
Распределение гидродинамических параметров и структура потоков в барботажном слое на ситчатой тарелке
В программной среде моделирования Компас создана модель ситчатой тарелки диаметром 700 мм и выполнено гидродинамическое моделирование в вычислительной среде Fluent системы Ansys.
Область барботажного слоя покрыта трехгранными конечными элементами, число которых приведено в таблице 1.
Для моделирования двухфазной среды используется Multiphase Model (Эйлера). В качестве сплошной фазы используется вода, в качестве дисперсной фазы - воздух. Задается шаг дискретизации по времени At = 10 "4
личество итераций на каждом шаге n Таблица 1 - Параметры сетки
с, число шагов и минимальное ко-0 = 50.
Количество трехгранных элементов в 18209087
области слоя
Количество трехгранных элементов 687788
на стенке
Количество трехгранных элементов 117380
на входе газа
Количество трехгранных элементов 43
на входе жидкости
Количество трехгранных элементов 2017
на выходе газа
Количество трехгранных элементов 41
на выходе жидкости
Количество трехгранных элементов 6693
на оси симметрии
Для инициализации решения вводятся начальные значения параметров (табл. 2).
Скорость газа на входе, м/с 5 8 8
Скорость жидкости на входе, м/с 0,043 0,06 0,044
Скорость диссипации энергии (1), м2/с3 2,99 3,18 3,2
Коэффициент турбулентной вязкости, м2/с (2) 0,001 0,00233 0,00232
Давление на входе газового потока, Па (3) 504 579 549
Давление на выходе газового потока, Па 0,33 0,33 0,33
Площадь межфазной поверхности, отнесенная к объему жидкости, 1/м (4) 287 150 150
Масштаб пульсаций, м (4) 0,0044 0,0813 0,008
Кинетическая энергия, м2/с2 (5) 0,183 0,287 0,287
Начальные условия получены в результате расчета параметров по приведенным выражениям [3 - 6]:
So^o
ы;
Рг +p*ghc
SkPr^3
£ =
и*ж =
Рг
РжЧк
(1)
vT = _И = 1.1и*ж ,(2) Рж s
\ 2Рж
P = AFcyx+AFp +АРпж, АРсух =С^2Р^2 , AFCT = 4ст/d0, АРпж = 1.3Ипж(1- Ф)Ржд ,
So
Aa
VЖ
( 2 Рг ю0
2
+Р жghc
^ SkPrю2 2
VЖ(pгU,г + 2 ст/ R3) k =
I = ■
\
и4 1.1^ с,,
C,= 0.009.
(3)
(4)
(5)
Здесь е - скорость диссипации энергии, So.Sk -свободное сечение тарелки и колонны, ко юк - скорость газа в отверстиях и в колонне, рг,рж - плотность газовой и жидкой фаз, и*ж - динамическая скорость трения жидкости на границе раздела фаз, hф - высота газового факела, hст - высота столба
жидкости на тарелке, V т - коэффициент турбулентной диффузии, ст - поверхностное натяжение, ф -газосодержание, £ - масштаб турбулентных пульсаций, к - кинетическая энергия среды, V■¡к - объем жидкости на тарелке, Аа - площадь поверхности переноса импульса.
По окончании расчета (число итераций на последнем шаге дискретизации п « 5 ) проводится вывод результатов в графическом или текстовом виде.
При экспериментальном исследовании структуры потоков на ситчатой тарелке диаметром 700 мм в работах Кафарова, Комиссарова, Амангалиева, Тё установлено наличие зон с различной степенью перемешивания, байпасирования и рециркуляции по-
2
v
т
U
ж
токов. За приемным порогом образуется циркуляционный контур, такой же контур с противоположным вращением образуется перед сливным порогом.
К таким же результатам приводит выполненное моделирование гидродинамических параметров методом вычислительной гидродинамики. На рис. (1 -3) показано распределение гидродинамических характеристик в барботажном слое на ситчатой тарелке по длине пути жидкости х, в поперечном направлении у и высоте слоя z.
На рис. 1 показано распределение продольной (вдоль пути жидкости) скорости и наличие трех характерных областей: за приемным порогом - область отрицательной скорости, обратного циркуляционного движения, у сливного порога область циркуляционного движения противоположного направления, и промежуточная область, ядро потока со слабо растущей продольной скоростью. Размеры этих областей, как и в экспериментах школы Кафа-рова составляют: ширина области перемешивания за приемным порогом и ширина области перемешивания у сливной планки около 20% от длины пути жидкости. При увеличении скорости газа в отверстиях тарелки с 5 до 8 м/с и расхода жидкости ширина зоны полного перемешивания у сливного порога уменьшается незначительно.
По высоте барботажного слоя (по координате z ) наблюдаем две характерные области: область ограниченная плоскостью тарелки (z = 0) и высотой факела, которая составляет около 10 мм, где расходуется основная часть кинетической энергии газа и наблюдаются высокие градиенты скорости жидкости и газа.
Продольная составляющая скорости их в этой области отрицательна во всех точках по х = 50,150, 350, 450 и у = 0, 260 (вдоль оси и цар-ги) (рис. 3). Следовательно, по дну тарелки от сливного порога к приемному рециклом возвращается часть жидкости. Расход этой жидкости определяется по скорости в середине сечения потока (hф /4).
Доля рецикла R определяется в виде отношения R = Lр / L, где Ьр - расход рецикла Lр = ирИфПр / 4,
где hф - высота факела, Пр - периметр потока рецикла, ир - средняя скорость рецикла.
их, м/с
Рис. 2 - Изменение продольной скорости по высоте барботажного слоя: 1 - за приемным порогом при х=50, у=260 (мм); 2 - у сливного порога при х=450, у=260 мм; 3 - в ядре потока в сечениях при х=150 (мм), у=260 (мм)
Uz , м/с
Рис. 1 - Изменение скорости в продольном направлении вдоль оси тарелки (у=0 мм) в различных сечениях по высоте слоя z (мм)
Рис. 3 - Изменение вертикальной составляющей скорости по высоте барботажного: 1 - за приемным и у сливного порогов при х=50; 450 (мм), у=0; 260 (мм); 2 - в ядре потока в сечениях при х=150;350 (мм), у=0; 260 (мм)
Скорость газа стремительно падает от входа на тарелку до распада газового факела на пузыри. По профилю скорости газа определяется высота факела Иф «8мм.
За пределами этой области (выше высоты факела) в промежуточной зоне и зонах полного перемешивания по высоте слоя z все три компоненты вектора скорости практически постоянны (рис. 3), что подтверждается экспериментальными данными о полном перемешивании жидкости по высоте барбо-тажного слоя. Продольная составляющая скорости ux постоянна в поперечном направлении практически до стенки и только в её окрестности сказывается влияние стенки.
Поперечная составляющая скорости Uy «^ и
близка к нулю. В этом случае можно считать близким к нулю конвективный поток массы uy 5x1 5у К 0 .
Изменение вертикальной составляющей скорости по высоте барботажного слоя 5uz/ 51 в области ядра потока, в силу полного перемешивания, равно нулю.
z, м
z, м
X, м
На основании проведенного анализа структуры потоков жидкости и газа на тарелке составляется комбинированная модель структуры потока (рис. 4).
Lx0 Зона (1 + R)Lx1 -Ядро (1 + R) Lx д Зона Lxebix
ПП-1 потока ПП-2
RLx вых
Рис. 4 - Структурная схема комбинированной модели структуры потока жидкости
На основе структурной схемы комбинированной модели составляются уравнения математической модели.
В зоне полного перемешивания ПП-1 уравнение материального баланса примет вид:
LXo + RLxвыx = (1 + R)Lx1. (6)
В ядре потока уравнение массопереноса представляется в виде:
dxn
ux^- = DT x d£, T
d2xr
d£
2
(7)
с граничными условиями:
Хд = х! при % = %ПП1, = 0 пРи % = %Д •
В зоне полного перемешивания III1-2 уравнение материального баланса запишем в виде:
(1 + Р)1_Хд = (1 + Р)1_Хвы,х (8)
Здесь % - продольная координата вдоль пути движения жидкости; %ПП1> %д - координата выхода из зоны ПП-1 и диффузионной зоны; хо, хвых, Х1, Хд -
концентрации жидкости на входе, выходе из тарелки, в области перемешивания ПП-1 и диффузионной зоне; D т « v т - коэффициент турбулентной диффузии.
Коэффициент турбулентной диффузии DT » vT, скорость жидкости в в ядре потока u x , границы зон £ ПП1, £д , доля рецикла R определяются по результатам гидродинамического моделирования.
Предлагаемый подход моделирования структуры потока жидкости и газа на контактных устройствах средствами вычислительной гидродинамики Ansys Fluent без привлечения экспериментальных исследований на физических макетах различного размера позволяет существенно сократить материальные затраты и сроки разработки контактного устройства.
Литература
1. Е.Н. Судаков, Расчеты основных процессов и аппаратов нефтепеработки. Химия, Москва, 1979. 568 с. 2. Ю.И. Дытнерский. Основные процессы и аппараты химической технологии. Альянс, Москва, 2007. 496с.
2. В.И. Елизаров, Д.В. Елизаров, С.А. Мерзляков, С.Г. Дьяконов, Теор. основы хим. технол., 46, 5, 483-490 (2012).
3. С.Г. Дьяконов, В.В. Елизаров, В.И. Елизаров Теоретические основы проектирования промышленных аппат-ратов химической технологии на базе сопряженного физического и математического моделирования. КГТУ, Казань 2009. 456 с. 3.
4. Д.В. Елизаров, В.И. Елизаров, С.А. Мерзляков, Вестник Казан. технол. ун-та, 16, 12, 206-210 (2013).
5. Д.В. Елизаров, В.В. Елизаров, С.А. Мерзляков, Вестник Казан. технол. ун-та, 17, 24, 300-303 (2014).
© С. А. Мерзляков - кандидат технических наук, доцент кафедры Автоматизации технологических процессов и производств, Нижнекамский химико-технологический институт (филиал) ФГБОУ ВО «Казанский национальный исследовательский технологический университет», Нижнекамск, Россия, sergew@inbox.ru; Д. В. Елизаров - кандидат технических наук, доцент кафедры Автоматизации технологических процессов и производств, Нижнекамский химико-технологический институт (филиал) ФГБОУ ВО «Казанский национальный исследовательский технологический университет», Нижнекамск, Россия; В. И. Елизаров - доктор технических наук, профессор кафедры Автоматизации технологических процессов и производств, Нижнекамский химико-технологический институт (филиал) ФГБОУ ВО «Казанский национальный исследовательский технологический университет», Нижнекамск, Россия; С. Г. Дьяконов - доктор технических наук, профессор, Советник при ректорате ФГБОУ ВО «Казанский национальный исследовательский технологический университет», Казань, Россия.
© S. Merzlyakov - Ph.D., Associate professor, Department of automation of technological processes and production, Nizhnekamsk chemical -technological Institute (branch) FGBOU VO "Kazan State technological University», Nizhnekamsk, Russia, sergew@inbox.ru, D. Elizarov - Ph.D., Associate Professor, Department of automation of technological processes and production, Nizhnekamsk Chemical-Technological Institute (branch) FGBOU VO "Kazan State Technological University», Nizhnekamsk, Russia, V. Elizarov - Doctor of Technical Sciences, professor of the Department of automation of technological processes and production, Nizhnekamsk chemical-technological Institute (branch) FGBOU VO "Kazan State technological University», Nizhnekamsk, Russia, S. Djakonov, Doctor of Technical Sciences, professor Advisor to the rector FGBOU VO "Kazan State Technological University", Kazan, Russia.
