Расчет аэродинамики циклонной камеры

2011

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Математика и механика

№ 1

МЕХАНИКА

УДК 531.351

М.А. Бубенчиков, И.А. Иванова

РАСЧЕТ АЭРОДИНАМИКИ ЦИКЛОННОИ КАМЕРЫ

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

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

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

Уместным будет сделать еще одно вводное замечание. В настоящей работе рассмотрены изотермические течения газа, реализующиеся при дозвуковых скоростях и перепадах давления порядка 10-2 атм. В этих условиях газ ведет себя как несжимаемая жидкость. Поэтому в дальнейшем будет использоваться термин «несжимаемая среда» вместо «газ» или «воздух».

Для случая ламинарного течения (Яе = Пср - Аг/V < 2000) несжимаемой ньютоновской жидкости с постоянной вязкостью в отсутствие массовых сил уравнения Навье - Стокса и уравнение неразрывности в цилиндрических координатах (г,

д д

ф, г) для двумерного стационарного осесимметричного случая (— = 0 , — = 0,

дґ дф

Ш Ф 0) течения могут быть представлены в форме [1]:

Система определяющих уравнений

дг дг р дг

дУ_ + у д¥_ _ ^ = _

дг дг г р дг

dW dW VW f 2 W'

и— + V— +----------= v[V2W-— I ;

dz dr r v r

U+1 d

dz r dr

,2 d2 d2 1 d

Здесь V =—- + —- + .

дг2 дг2 г дг

Уравнения (1) - (4) являются эллиптическими, поэтому для их однозначного разрешения по краю ограниченного куска плоскости, на котором они определены, необходимо выставить граничные условия для и, V и Ж. В общем случае поставленная таким образом задача может быть решена только численно. В динамике вязкой жидкости известна проблема расчета распределения давления, согласованного с полем скорости. Чтобы избежать этой проблемы, в исходной системе от естественных переменных (и, V, Ж, р) перейдем к переменным «функция тока - завихренность» (у, 4, Ж):

1 д^ =и 1 д^=^ ?=—-— (5) г дг , г дг , дг дг

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

=

V г дг )

и^Ъ-И-Ж^ =М-4);

дг дг г г дг V г 2)

дЖ дЖ VW ( 2 Ж \

и — + V— + — = v\v2W---\.

дг дг г V г

К этим уравнениям следует добавить соотношения (5).

Физическая область течения

Так как течение является осесимметричным и независимых переменных только две - z и r, то можно рассматривать любое осевое сечение цилиндрической камеры, а область интегрирования в этом случае будет плоской фигурой, в нашем примере прямоугольником.

Граничные условия

Все представленные ниже граничные условия для искомых величин получены из физических условий для компонент скорости:

UL = Ulr„ = 0, Vlr„ = Увх = const, ^ = Уы = cons,.

г 4

Рис. 1. Продольное сечение циклонной камеры. Сплошные линии представляют непроницаемые границы, пунктирные - проницаемые

Граничные условия для функции тока у следующие:

Мг. = 0 ( ^ < = бвх , Чги = у—  бвх , Чгвых

у

бвых .

Здесь (вх = -Увх г1Квх , Ьвх - протяженность входной зоны, Увх - величина скорости жидкости на входе (постоянная величина), (вых = -Увыхг1^вых, ¥вых - величина скорости жидкости на выходе (постоянная положительная величина), Увых - протяженность выходной зоны.

Так как имеется лишь один вход и один выход, то условие интегрального баланса массы будет заключаться в выполнении условия (вых = (вх.

Граничные условия для завихренности £:

^1 =

дУ_

дг

ди

дг

д¥

дг

г г г 1

ди

дг

Граничные условия для окружной скорости:

Щ г1 =аг (. ^ Щ гвх =0,

дШ

г

= 0.

Задание численных граничных условий оказывает существенное влияние не только на устойчивость, но и на точность решения конечного разностного уравнения для завихренности и задачи в целом. Поскольку и и V не определяются непосредственно из системы (6) - (8), а находятся после ее решения как разностные производные от распределений у, то и в разностном аналоге граничных условий для £ на стенке следует перейти от численных производных от компонент скорости (см. (12)) к соответствующим производным от функции тока. Имеется множество вариантов записи таких выражений, однако простейший выглядит следующим образом:

2

='

Ап

- + о (Ап).

Здесь, независимо от ориентации стенки, уш - значение функции тока на самой стенке, уш+1 - в первой узловой точке, расположенной в жидкости, Ал - расстояние по нормали от стенки до ближайшей узловой точки.

г

Такое условие первого порядка было впервые предложено в 1928 г. в работе Тома [2] и широко используется до настоящего времени. Это условие очень надежно и часто приводит к результатам, достаточно хорошо согласующимся с результатами, полученными при помощи форм высших порядков граничного условия для вихря [3].

Разностный вид граничных условий для завихренности £ будет следующим:

ЄІГ =

д¥

дг

2 у(2, } )-У(1,1) •О) Аг2

+о(Аг).

ди

2 у

+о(Аг),

г=і Г

= - ди

г

Аг

+о(Аг),

2 у(і,2)-у(і,1) '(1) Аг2

+ о(Аг).

Здесь К - количество интервалов разбиения области интегрирования в г-направлении, N - в г-направлении.

Метод решения

Решение уравнений (6) - (8) с граничными условиями (11) - (13) будем строить с использованием конечных разностей и итерационного метода последовательной перерелаксации (ППР) [4]. При аппроксимации диффузионных членов будут применены симметричные разности. Аппроксимация конвективных членов будет проведена с использованием разностей против потока. Дифференциальные члены уравнения Пуассона для у заменяются симметричными разностями, имеющими второй порядок аппроксимации:

Л ^-2у,., +4,+1, .

^ д2^ + г д ( 1 ду

аг 2

дг V г дг

Аг 2

+г.

а, +1/2 (+1 -Уі,і )- а,-1/2 (( -Уг,}-1) Аг 2

+ о

(2, А).

Здесь а

І-1/2

г-1 + О І-1

-7+1/2

г + г+1 І+1

Аппроксимация конвективных членов проведена следующим образом:

и-

,дП_

дг

тт+ иі, І иі-1,І т т — иі+1,І иі, І /а Ч

= и+і— -------------— + иі і------- -------— + о (Аг).

і, і

Аг

і і

Аг

где иг,] =

ии + Р,А

ии =

ии - и,А

2 2 Аналогично для конвективного переноса в г-направлении:

V

ди

дг

тл+ иі,І иі,І-1 т /— иі, І+1 иі, І /  \

= у+— ----------- — + V ,— -----------------— + о (Аг).

і, і

Аг

Аг

V  + V . V . -V \

О Т7+ '• 1 '• Л 17- '•1 '• Л

Здесь ¥1} =------^, ¥1} =----------

Приравнивая разностное выражение (14) значению -^ц и выражая из полученного равенства у, найдем рекуррентную формулу для пересчета сеточных значений функции тока:

{

^ __ У дг,./ ./_____т г 1./_____т ^1 У /_У \ ] 1/2 т г,./ 1________________уі-1/2 т г, уі-1 /

где введено обозначение т = .

Дг 2

Рекуррентные формулы для пересчета значений £у и Ж,-,,- будут выглядеть следующим образом:

^',1 = Г-Дг2Р

Здесь Р _У

( г. г. , \

2 2 1

Аг2 Аг 2 гі 2 .

V У У

+И.+1^1

Аг Аг

е _ + ЦIу -1,у _ и;, Ъ +1,у + У*-1 _ ^+1 + «у Г «+1,у - «-Ц

Г/ Аг Аг Аг Аг г}^ Аг

-1/2»;,у-1 + Г,.„#№1)+)+у.

где а У » ! + Ц+»;-1,у - иУ-ы + У'«»-1 - +1,

г- Аг Аг Аг Аг

г _ г--1 + г г _ г + гу+1

0-1/2 - , , г+1/2

22

Последовательное применение формулы (15) в движении по внутренним сеточным узлам предполагает использование новых значений у там, где процесс перебора узлов уже завершен, и значений у на прошлом итерационном слое для узлов, еще неподверженных пересчету. Вслед за тем как искомая величина последовательно перевычислена во всех внутренних узлах, она перевычисляется в граничных узлах. После того, как все искомые величины таким образом найдены, одна глобальная итерация считается законченной. Выполнив порядка 103^105 глобальных итераций, мы завершаем расчет.

Результаты расчетов

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

Во-первых, это тривиальный тест: QвЬIX = 0вх = 0 , п = 0, где п - частота вращения барабанов. В результате рассчитываются нулевые распределения всех величин.

Во-вторых, тест твердотельного вращения газа в камере: и = V = 0, п Ф 0, имеем линейное распределение окружной скорости по толщине вращающегося слоя.

Третий тест - прямоточное течение среды в кольцевом канале. Получаем автомодельный профиль продольной компоненты скорости [5].

Четвертый тест: QвЬIX = Qвх = 0, и Ф 0, V Ф 0, Ж Ф 0 (п Ф 0). Генерируются

вихри Тейлора.

На рис. 2. приведены примеры расчетов течения газа, когда входное и выходное кольцо равны по размерам, и соответственно скорости радиальной подачи и выхода газа тоже равны: Ивх = Vвых = 1,0 м/с . Примеры отличаются интенсивностью крутки камеры. Как видим из рисунка, с ростом п характер течения становится более сложным, появляется больше тороидальных вихрей с более высокой интенсивностью движения частиц газа в этих образованиях. Однако с ростом п

0,068

0,066

0,064

0,062

0,060

г~— —

0,1

0,2

а

0,3

0,068

0,066

0,064

0,062

0,060

тЙШ

0,1

0,2

С

0,3

0,068 0,066 -0,064 0,062 0,060

\

/И

ш>м

т,/ Ш \ , , «

^ \&\\ ^ \ II )

Ж

0,1

0,2

Ь

0,3

0,068  0,066  0,064  0,062 0,060

0,1

0,2

а

0,3

Рис. 2. Линии тока, полученные при а - п = 10об/с; Ь - п = 50об/с;

различных величинах крутки камеры: с - п = 100об/с; а - п = 200об/с

г, м

г, м

0

0

г, м

г, м

г, м

г, м

0

0

г, м

г, м

увеличивается также количество глобальных итераций K, обеспечивающих сходимость вычислительного процесса. Так, при n = 10 K* = 103; n = 50 K* = 104; n = 200 K* = 105.

Заключение

Как показали вычисления, предложенный алгоритм сходится даже при числах Рейнольдса, формально превышающих критическое число перехода к турбулентному режиму, и позволяет проводить расчеты при частотах вращения барабанов, отвечающих режиму выделения наночастиц из воздуха, т.е. при n ~ 100 об/с и более. Однако количество глобальных итераций, обеспечивающих сходимость, существенно зависят от частоты вращения барабанов. Таким образом, в настоящей работе предлагается простой способ вычислений, позволяющий эффективно решать задачи о закрученных течениях вязкой жидкости без привлечения сверхсложных вычислительных технологий, использования дорогостоящих программных пакетов и больших вычислительных мощностей ЭВМ.

ЛИТЕРАТУРА

1. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: учеб. пособие: в 10 т. М.: Наука, 1986. Т. VI: Гидродинамика. 736 с.

2. Thom A. An investigation of fluid flow in two dimensions // Aerospace Research Center, К and M. 1928. No. 1194.

3. Роуч П. Вычислительная гидродинамика. М.: Мир, 1983. 648 с.

4. Численные методы в динамике жидкостей / Э. Джеймсон [и др.]. М.: Мир, 1981. 408 с.

5. Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1973. 607 с.

Статья принята в печать 04.02.2011г.

Bubenchikov M.A., Ivanova I.A. CALCULATION OF AERODYNAMICS VORTEX CHAMBER. In this paper, the simple iteration method usually used to solve diffusion problems is generalized to the solutions of transport equations with convection. The developed computational technology is used to analyze the effect of camera rotation intensity on the nature of the movement of the gaseous medium against the walls of the rotating section of the apparatus.

Keywords: gas centrifuging, vortex chamber, mathematical modeling, viscosity, stream function, vorticity, method of simple iteration.

BUBENCHIKOV Mikhail Alekseevich

E-mail: michael121@mail.ru

IVANOVA Irina Aleksandrovna

E-mail: vestnik_tgu_mm@math.tsu.ru