Scientific journal
Modern problems of science and education
ISSN 2070-7428
"Перечень" ВАК
ИФ РИНЦ = 1,006

SIMULATION OF TURBULENT GAS-DISPERSION FLOWS

Sugak E.V. 1 Sugak A.V. 2
1 Siberian State Aerospace University
2 Yaroslavl State Technical University
A new approach to the modeling of turbulent gas-dispersed flows is considered. The technique of probabilistic and statistical modeling of gas-dispersed flows incorporating the deterministic-stochastic nature of the turbulent motion of the continuous and dispersed phases with use of probability theory and stochastic processes. A model of turbulent gas-dispersion flow in a cylindrical channel, the method of calculation of the concentration profile of the particles and changes along the channel. The developed methods for simulation of turbulent flows of gas-dispersion can be used in the analysis, the intensification and improvement of the hydrodynamic, heat and mass transfer processes in heterophase system is in intensive interaction phases of modeling and calculation processes and devices of chemical technologies for purification of industrial waste gases from gaseous impurities and highly apparatuses with intense hydrodynamic regimes and dust filters.
turbulence
processes and devices
aerosols
two-phase systems
aerodynamics

Во многих энергетических, тепло- и массообменных процессах и аппаратах используются режимы интенсивного взаимодействия веществ в различных агрегатных состояниях. Интенсификация таких процессов часто связана с повышением скорости движения фаз, переходом в существенно турбулентные режимы течения [1–4,7,9,10].

Течение двухфазных турбулентных газодисперсных потоков в реальных аппаратах обычно настолько сложно, что построение их полного и точного математического описания без каких-либо упрощающих допущений не представляется возможным [1,3,10]. При этом возникают трудности как физического, так и математического характера, связанные как с высокой размерностью задач, так и с целым спектром эффектов взаимодействия фаз разной природы, значительная часть которых носит вероятностно-статистический характер. В связи с этим методы теоретических и экспериментальных исследований аэродисперсной динамики развиты в гораздо меньшей степени, чем методы классической гидроаэродинамики [1,2].

Для описания гидродинамики двухфазных потоков существует два основных подхода - континуальный и статистический [1–3,10]. Несмотря на существенные различия, основной общей проблемой для них является необходимость учета большого количества разнообразных факторов и процессов межфазного взаимодействия, имеющих различную физическую природу. Характер движения дисперсной фазы в турбулентном потоке носит вероятностно-статистический характер, и попытки его описания детерминированными зависимостями существенно снижают возможности анализа и оптимизации технологических процессов и оборудования. Использование детерминированных методов в большинстве случаев позволяет определять только ориентировочные или усредненные значения параметров и характеристик процесса, что часто приводит к ошибкам, снижению точности расчетов или необходимости введения эмпирических коэффициентов, не имеющих четкого физического смысла [4,10]. Недостатки существующих моделей движения дисперсной фазы в газодисперсных потоках в трубах и каналах могут быть устранены с использованием вероятностно-статистического подхода к их моделированию с учетом некоторых упрощающих допущений [4–8].

Если траекторию движения частицы в газовом потоке рассматривать как суммарный случайный путь, то в любой момент времени любая из координат частицы может быть представлена как сумма детерминированной и случайной составляющих. Детерминированная составляющая может быть описана законами классической механики и механики сплошной среды, случайная составляющая описывается методами теории случайных процессов [5–8].

Если рассматривать движение частицы под действием турбулентных пульсаций среды как последовательность скачкообразных перемещений длиной l через малые промежутки времени Δt в одном из шести возможных направлений в ортогональной системе координат x-y-z, то траектория движения будет представлять собой трехмерную ломаную, а направление движения в каждый момент времени будет определяться соответствующими вероятностями pi: p+x, p-x, p+y, p-y, p+z и p-z, при этом в любой момент времени p+x(t) + p-x(t) + p+y(t) + p-y(t) + p+z(t) + p-z(t) = 1. В отсутствие конвективного движения и влияния внешних сил (или в системе координат, движущейся с газом) частица совершает только случайные движения и при изотропной турбулентности все направления движения равновероятны и pi = 1/(2s) (где s - размерность системы координат, для трехмерной системы s = 3 и pi = 1/6).

При течении в канале можно перейти к двухмерной системе координат (s = 2) и рассматривать движение частиц в одном из четырех возможных направлений - два вдоль поверхности канала (+z и -z) и два в поперечном направлении (+y и -y). Тогда траектория частицы будет представлять двухмерную ломаную, а направление движения в каждый момент будет определяться вероятностями p+z, p-z, p+y и p-y, при p+z > p-z = p+y= p-y и p+z+p-z+p+у+p-y = 1 (рис. 1).

Рис. 1. Схема и пример расчета движения частицы

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

Рис. 2. Дискретное вероятностное моделирование траекторий 10 частиц

Однако можно заметить, что при таком подходе в любой момент времени частица может находиться в одном из узлов сетки с ячейками со сторонами lz´ly. Каждое из таких положений можно рассматривать как возможное состояние частицы в момент времени t, которое характеризуется соответствующей вероятностью P(i,j,t) (где i и j - номера узлов решетки, т.е. безразмерные координаты частицы: i = yi/ly, j = zi/lz).

Пусть известны вероятности всех положений частицы в момент времени t. Рассмотрим изменение вероятности нахождения частицы в положении (i,j)) через малый промежуток времени Δt. При этом интервал Δt предполагается настолько малым, что вероятность того, что в течение этого времени произойдет более одного перехода пренебрежимо мала.

В момент времени t + Δt частица может находиться в положении (i,j) только в одном из двух случаев (рис. 3):

1. В момент t частица находилась в одном из четырех соседних положений (i,j-1), (i,j+1), (i-1,j), (i+1,j) и в течение времени Δt произошел один из переходов: (i,j-1)→(i,j), (i,j+1)→(i,j), (i-1,j)→(i,j) или (i+1,j)→(i,j). Суммарная вероятность этих событий по теореме сложения вероятностей

P1(i,j,t+Δt) = p+z(i,j-1,Δt)P(i,j-1,t) + p-z(i,j+1,Δt)P(i,j+1,t) + p+y(i-1,j,Δt)P(i-1,j,t) + p-y(i+1,j,Δt)P(i+1,j,t), (1)

где P(i,j-1,t), P(i,j+1,t), P(i-1,j,t) и P(i+1,j,t) - вероятности нахождения частицы в момент времени t в положениях (i,j-1), (i,j+1), (i-1,j) и (i+1,j), соответственно; p+z(i,j-1,Δt), p-z(i,j+1,Δt), p+y(i-1,j,Δt) и p-y(i+1,j,Δt) - вероятности переходов в положение (i,j) из положений (i,j-1), (i,j+1), (i-1,j) и (i+1,j), соответственно.

2. В момент времени t частица уже находилась в положении (i,j) и в течение времени Δt не перешла ни в одно из соседних положений, т.е. не состоялись переходы (i,j)→(i,j-1), (i,j)→(i,j+1), (i,j)→(i-1,j) и (i,j)→(i+1,j). Вероятность такого события по теореме умножения вероятностей

P2(i,j,t+Δt) = [1 - p-z(i,j,Δt)][1 - p+z(i,j,Δt)][1 - p+y(i,j,Δt)][1 - p-y(i,j,Δt)]P(i,j,t). (2)

где p+z(i,j,Δt), p-z(i,j,Δt), p+y(i,j,Δt) и p-y(i,j,Δt) - вероятности переходов из положения (i,j) в соответствующие соседние положения.

 

Рис.3. Фрагмент графа переходов частицы в потоке

Тогда полная вероятность того, что частица в момент времени t+Δt будет находиться в положении (i,j) по теореме сложения вероятностей

P(i,j,t+Δt) = P1(i,j,t) + P2(i,j,t) = p+z(i,j-1,Δt)P(i,j-1,t) + p-z(i,j+1,Δt)P(i,j+1,t) + p+y(i-1,j,Δt)P(i-1,j,t) + p-y(i+1,j,Δt)P(i+1,j,t) +

+ [1 - p-z(i,j,Δt)][1 - p+z(i,j,Δt)][1 - p-y(i,j,Δt)][1 - p+y(i,j,Δt)]P(i,j,t). (3)

Для рекуррентного пуассоновского потока событий при li(i,j)Δt << 1 [8]

pi(i,j,Δt) = 1 - exp[-li(i,j)Δt] » 1 - [1 - li(i,j)Δt] = li(i,j)Δt, (4)

где li (i = +z,-z,+y,-y) - интенсивности соответствующих переходов, которые в турбулентном потоке определяются частотой турбулентных пульсаций газа, с-1.

Аналогичные приближения можно записать для всех соседних положений. Тогда, если пренебречь величинами второго порядка малости, уравнение (3) примет вид

P(i,j,t+Δt) » [l+z(i,j-1)P(i,j-1,t) + l-z(i,j+1)P(i,j+1,t) + l+y(i-1,j)P(i-1,j,t) +

+ l-y(i+1,j)P(i+1,j,t)]Δt + {1 - [l+z(i,j) + l-z(i,j) + l+y(i,j) + l-y(i,j)]Δt}P(i,j,t). (5)

Из выражения (5) можно получить

[P(i,j,t+Δt) - P(i,j,t)]/Δt = l+z(i,j-1)P(i,j-1,t) + l-z(i,j+1)P(i,j+1,t) + l+y(i-1,j)P(i-1,j,t) +

+l-y(i+1,j)P(i+1,j,t) - [l+z(i,j) + l-z(i,j) + l+y(i,j) + l-y(i,j)]P(i,j,t). (6)

Переходя к пределу при Δt→0, получим дифференциальное уравнение относительно функции P(i,j,t):

l+z(i,j-1)P(i,j-1,t) + l-z(i,j+1)P(i,j+1,t) + l+y(i-1,j)P(i-1,j,t) +

+ l-y(i+1,j)P(i+1,j,t) - [l+z(i,j) + l-z(i,j) + l+y(i,j) + l-y(i,j)]P(i,j,t). (7)

Если известны значения интенсивностей переходов li (i = +z, -z, +y, -y), то, составив для каждого из возможных положений частицы уравнение вида (7), получим систему дифференциальных уравнений. Так как полученная система линейно зависима, то одно из уравнений необходимо заменить нормирующим условием вида

. (8)

Вероятность нахождения одиночной частицы P(i,j,t) в любом из положений в соответствии с законом больших чисел означает долю частиц N(i,j,t)/N от их общего числа в системе, находящихся в объеме V(i,j) сечением lz´ly в момент времени t, то есть

P(i,j,t) = N(i,j,t)/N = n(i,j,t)V(i,j)/N, (9)

где n (i,j,t) - локальная численная концентрация частиц, м-3.

Тогда уравнение (7) принимает вид

l+z(i,j-1)n(i,j-1,t)V(i,j-1) + l-z(i,j+1)n(i,j+1,t)V(i,j+1) + l+y(i-1,j)n(i-1,j,t)V(i-1,j) + l-y(i+1,j)n(i+1,j,t)V(i+1,j) -

- [l+z(i,j) + l-z(i,j) + l+y(i,j) + l-y(i,j)]n(i,j,t)V(i,j). (10)

При постоянном расходе и установившемся течении характеристики турбулентных пульсаций газа (амплитуду и частоту) можно считать постоянными в каждой точке и одинаковыми вдоль направления течения. Тогда V(i,j-1) = V(i,j+1) = V(i,j) = V(i), V(i-1,j) = V(i-1), V(i+1,j) = V(i+1) и в уравнении (10)

, (11)

и оно принимает вид

l+z(i)[n(i,j-1,t) - n(i,j,t)] + l-z(i)[n(i,j+1,t) - n(i,j,t)] +

. (12)

Интенсивности переходов li определяются детерминированными скоростями и интенсивностями турбулентных пульсаций в соответствующих направлениях. Скорость частицы в любом направлении определяется суммой детерминированной и случайной составляющих: ux =`ux + ux¢. Тогда общая интенсивность переходов вдоль любой из осей определяется суммой интенсивности переходов с детерминированной осредненной скоростью вдоль этой оси`ux/lx и интенсивности турбулентных пульсаций вдоль той же оси lx: l±x =`ux/lx ± lx.

Если считать, что при течении аэродисперсной смеси в канале в отсутствие внешних сил детерминированная осредненная скорость частиц направлена строго вдоль поверхности по оси z, то в уравнении (12) l+z(i) =`u(i)/lz(i) + lz(i), l-z(i) = lz(i), l+y(i) = l-y(i) = ly(i), и оно примет вид (знак усреднения продольной осредненной скорости частиц далее опущен)

u(i)[n(i,j-1,t) - n(i,j,t)]/lz(i) - 2(lz + ly)n(i,j,t) + lz(i)[n(i,j-1,t) +n(i,j+1,t)] +

+ ly(i)[n(i-1,j,t)ly(i-1)lz(i-1) + n(i+1,j,t)ly(i+1)lz(i+1)]/[ly(i)lz(i)]. (13)

Система дифференциальных уравнений вида (13) при заданных начальных и граничных условиях дает возможность определить текущую локальную концентрацию частиц и ее изменение. Однако в большинстве практически значимых случаев интерес представляет установившееся течение, при котором все параметры двухфазного потока постоянны. При этом ¶n(i,j,t)/¶t = 0, n(i,j,t) = n(i,j), и уравнение (13) принимает вид

. (14)

Система линейных алгебраических уравнений вида (14) дает возможность определить локальную концентрацию частиц при установившемся течении. Уравнения для узлов, расположенных на границах области (в начальном сечении и на стенке), в соответствии с заданными граничными условиями будут несколько отличаться от уравнения общего вида (14). Так на входе можно принять равномерное распределение частиц по сечению (n(i,0) = n0), а интенсивность прямых переходов в направлении течения) считать зависящей только от детерминированной скорости (p+z(i,0,Δt) = 1, l+z(i,0) = u(i)/lz(i)).

В канале с непоглощающими твердыми стенками можно считать, что после контакта со стенкой частицы остаются в непосредственной близости от нее. Это означает, что в любом граничном положении (i = 0, yi = 0) исключены как прямые, так и обратные радиальные переходы. Тогда для пристенных узлов уравнение (14) приобретает вид

. (15)

Совокупность уравнений (14) и (15) с учетом начальных и граничных условий образует систему, решение которой позволяет определить профиль концентрации частиц в любом сечении канала и его изменение по длине. Однако, учитывая масштаб турбулентных пульсаций, размерность системы будет очень большой, и расчет становится очень трудоемким.

Однако можно заметить, что расчетная схема на рис. 3 содержит точки, принадлежащие трем последовательным сечениям (слоям) по длине канала - j-1, j и j+1. Следовательно, если известны концентрации частиц двух предыдущих сечений, то можно рассчитать концентрацию частиц в следующем сечении. Соответственно, уравнения (14) и (15) могут быть представлены как итерационные уравнения вида

, (16)

, i = 1,2,…,k-1, j = 1,2,…,m. (17)

Турбулентную миграцию частиц можно интерпретировать как цепочку свободных инерционных пробегов с длиной ly(y) и lz(y), периодом T и частотой w = 1/T [4-6]. В уравнениях (16) и (17) можно также принять ly(i) = lz(i) = l(i), и они принимают вид

, j = 1,2,…,m, (18)

, (19)

Однако расчет по уравнениям (19) может быть затруднен из-за переменной по сечению канала длины инерционного пробега частиц. Так как, согласно принятой схеме, расстояние по осям между соседними узлами равно, соответственно, lz(i) и ly(i), то пристенные узлы будут «запаздывать» по сравнению с приосевыми, что приведет к искривлению расчетной сетки и деформированию формы ячеек. Поэтому целесообразно ввести сетку с постоянным по длине расстоянием между узлами Dz, равным, например, длине пробега частиц на оси канала. Тогда, предполагая функцию n(z) монотонной, при Dz << L, можно принять n(i,z+Dz) » n(i,z) + Dnz(i,j)Dz/lz(i). Для сокращения объема вычислений можно также принять допущение об изотропной турбулентности в каждой точке потока lz(i) = ly(i) = l(i) и lz(i) = ly(i) = l(i)). Тогда уравнение (19) примет вид [4]:

, (20)

или, переходя в обозначении узлов к нижним индексам, при допущении, что интенсивность турбулентных пульсаций газа и интенсивности переходов одинаковы по сечению канала (l(i) = li = l), то окончательно получим

, (21)

В случае осесимметричного потока для осевых узлов (при i = s) ls+1 = ls-1, ns+1 = ns-1 и уравнение (21) приобретает вид

. (22)

Уравнения (21) и (22) дают возможность последовательного расчета (от сечения к сечению) численной концентрации частиц и ее распределения по длине и сечению канала, исходя из интенсивности турбулентных пульсаций l и профилей осредненной скорости u(r) и длины пробега частиц l(r), которые могут быть получены на основании расчета течения, исходя из параметров потока и характеристик двухфазной среды [4]. Сочетание аналитического (для каждого микропроцесса) и численного (для процесса в целом) подходов позволяет упростить расчет (по сравнению с численным) и в то же время провести его более полный и точный анализ (по сравнению с аналитическим расчетом без учета случайных факторов).

Расчеты по приведенной методике показывают, что профиль концентрации частиц в турбулентном потоке характеризуется значительной неоднородностью, особенно в пристенной области (рис.4). Наиболее интенсивная перестройка профиля происходит на начальном участке канала, на котором профиль концентрации изменяется от выпукло-вогнутого до вогнутого (чашевидного), после чего меняется медленно, т.е. наступает динамическое равновесие, при котором число частиц, поступающих из приосевой зоны в пристенную, равно числу возвращающихся частиц. Длина участка стабилизации зависит, прежде всего, от размера частиц (от степени их увлечения турбулентными пульсациями) и при обычных условиях лежит в интервале от нескольких диаметров (для частиц диаметром менее 1 мкм - рис. 4а) до нескольких сотен диаметров (для частиц диаметром более 100 мкм - рис.4б). Установившийся профиль концентрации частиц при z→¥ можно определить из уравнений (21)-(22) при n(z+Dz) = n(z), из которых следует, что равновесная концентрация обратно пропорциональна квадрату длины инерционного пробега частиц (niр ~ li-2) [4].

а)

б)

Рис.4. Изменение профиля концентрации частиц по длине канала:

а -d = 1 мкм, б - d = 100 мкм, D = 51 мм, w = 10 м/с, r = 1000 кг/м3

Наличие различных по форме профилей концентрации дисперсной фазы в турбулентном потоке подтверждается экспериментальными исследованиями. В зависимости от степени увлечения профиль концентрации частиц в различных сечениях по длине канала может принимать разный вид - логарифмический, прямоугольный, чашевидный, седловидный или куполообразный [1,2,4,10].

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

Рецензенты:

Хакимзянов Гаяз Салимович, доктор физико-математических наук, профессор, ведущий научный сотрудник, Институт вычислительных технологий СО РАН, г. Новосибирск.

Чекалов Лев Валентинович, доктор технических наук, генеральный директор ЗАО «Кондор-Эко», Ярославская область, г. Семибратово.