Содержание

Введение                                                                                                                                         3

§ 1.     Физическая постановка задачи                                                                                        4

§ 2.     Математическая постановка задачи                                                                               7

§ 3.     Задача в безразмерном виде                                                                                           9

§ 4.     Алгоритм и его численная реализация                                                                          11

§ 5.     Результаты  расчетов                                                                                                      15

§ 6.     Сравнение с результатами двумерной модели                                                            21

Заключение                                                                                                                                 27

Список литературы                                                                                                                     27

Введение

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

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

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

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

Предметом исследования в настоящей работе является вихревая камера плазматрона. Целью работы является математическое моделирование течения газа в такой камере в трехмерном приближении. Ранее в [1] аналогичная задача решена в двумерном случае. Тем самым в данной работе предполагается уточнить результаты [1] и найти пределы применимости двумерной модели.

§ 1. Физическая постановка задачи

1.1           Физическая постановка задачи

На рис. 1 изображена схема рассматриваемой вихревой камеры.

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

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

1.2           Основные предположения:

1)    отверстия, через которые ведется подача газа, располагаются тангенциально и  симметрично относительно оси вихревой камеры;

2)    задача будет решаться для шести отверстий;

3)    в начальный момент времени газ покоится;

4)    на границе вдува задано постоянное давление , одинаковое для всех форсунок;

5)    через боковые отверстия газ поступает в основную часть форсунки с постоянной скоростью;

6)    расчеты будут вестись для дозвуковых течений  (числа Маха меньше 1);

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

8)    описание будет вестись только в эйлеровых переменных.

1.3           Обозначения

Будем использовать следующие обозначения: t – время;  - точка пространства;  - плотность газа;  - скорость газа в декартовой системе координат; - скорость газа в цилиндрических координатах ;  - объемная плотность полной энергии газа, где  - массовая плотность удельной внутренней энергии газа;  - давление;  - температура газа;  - вектор вихря;  - показатель адиабаты;  - универсальная газовая постоянная;  - область;  - граница области ;  – радиус вихревой камеры;  - радиус отверстия, через которое происходит вдув газа;  - скорость звука на входе; M - число Маха на входе;  - единичный тензор;  - тензор напряжений;  - коэффициент динамической вязкости;  - число Рейнольдса;  - начальная плотность;  - плотность на границе вдува;  - граничная скорость (скорость втекания газа);  - начальная температура (она же температура стенок); - давление на входе.

1.4           Критерий выбора оптимальной геометрии

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

                ,                        

                 ,                        

где  - единичные векторы, направленные вдоль осей .

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

§ 2. Математическая постановка задачи

2.1    Описание математической модели

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

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

Неизвестными величинами являются  - плотность,  = (u, v, w) – скорость газа,  - давление, T – температура,  - норма вектора вихря.

2.2                  Система уравнений газовой динамики

Ниже представлена система определяющих уравнений газовой динамики.

Система состоит из уравнений неразрывности, импульса и энергии.  Она дополнена уравнениями состояния для совершенного газа , начальными и граничными условиями (см. ниже). Задача решается в трехмерной «секторной» области (рис. 2).

Компоненты тензора напряжений  возьмем в виде [2, 4]:

                                   ,

где  -  коэффициент динамической вязкости.

Будем считать, что в начальный момент времени в области D выполнены следующие условия:

, , .

Ниже представлены основные граничные условия:

o       условие втекания: , , ,  в отверстии вдува ;

o       условие прилипания:  на твердых стенках камеры (в нашем случае это боковая часть, область около отверстия вдува и нижнее основание камеры), то есть на ;

o       условие свободного вытекания на верхнем основании камеры  (подробнее см. ниже);

o       условие периодичности на боковых границах сектора:

где  - угол сектора (в нашем случае ),  - угол между нижней гранью сектора и осью OX,  (в данном случае ) - нормальная к боковым граням составляющая скорости,  - радиальная составляющая скорости.

§ 3. Задача в безразмерном виде

3.1 Обозначения, математическая постановка и параметры задачи

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

Введем следующие безразмерные переменные: - безразмерная i–ая координата, отвечающая декартовой системе координат, - безразмерная i–ая координата скорости, , ,  - безразмерные время, плотность, давление, температура, внутренняя энергия и объемная плотность полной энергии соответственно,  - число Маха на входе, где - скорость звука на входе,  - число Рейнольдса. Тогда

где  - безразмерный тензор напряжений.

Начальные условия:

, , .

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

o       условие втекания: , , ;

o       условие непротекания: ;

o       условие свободного вытекания (см. ниже);

o       условие периодичности на боковых границах.

В размерных переменных радиус камеры  мм, начальная плотность  кг/м³, начальное давление  атм, начальная температура .   Скорость звука  м/c.

Расчеты проведены при разных числах Маха (0.2, 0.4, 0.6, 0.8, 1.0), и разных относительных радиусах отверстий вдува  (0.03, 0.04, 0.05, 0.06, 0.065). Для всех вариантов геометрических параметров (рис. 2) фиксированными остаются высота «сектора» , длина отверстия вдува равна .

Безразмерная начальная плотность газа . Безразмерная начальная температура газа . Безразмерная плотность газа на вдуве . Давление на входе в четыре раза больше начального давления . Показатель адиабаты в вычислениях взят равным . Число Рейнольдса Re = 548167.

Далее волны у безразмерных величин опущены.

Число Куранта для гиперболической части решаемой системы уравнений:

                                               ,

где  - максимальное собственное значение,  - шаг по времени,  - шаг по пространству. Для задачи в безразмерном виде  . Тогда справедлива оценка:

Число Куранта для параболической части решаемой системы уравнений:

                                               .

Для него справедлива оценка:

§ 4. Алгоритм и его численная реализация

4.1                        Программный комплекс

Разработан программный комплекс для решения задач трехмерной неидеальной газовой динамики в областях сложной геометрии, адаптированный под данную физическую задачу. Он позволяет получать решения, зависящие от безразмерных критериев, и помимо эффектов, свойственных идеальной газовой динамике, учитывать влияние вязкости. В основе алгоритма лежит прием «расщепления по физическим процессам» [10]. Численные решения, отвечающие различным физическим процессам, реализуются в отдельных программных модулях. Уравнения идеальной газовой динамики численно решаются методами, основанными на приближенном решении задачи Римана. Потоки через границы ячеек, отвечающие невязким и нетеплопроводным слагаемым, вычисляются методом Хартена-Лакса-ван Лира (HLLC) (см. [6, 7]). Алгоритм аппроксимации слагаемых, отвечающих эффектам вязкости и теплопроводности газа, основан на аппроксимации частных производных. Его идея взята из [8] и описана подробнее в [1]. На границе, на которой происходит  сводное вытекание, поставлено неотражающее граничное условие.

4.2                        Постановка неотражающих граничных условий

Более подробно постановка неотражающих граничных условий для систем гиперболических уравнений описана в [6, 12].

В нашем случае оператор решаемой системы уравнений представляется в виде суммы двух операторов – оператора гиперболической части системы и «параболического» оператора, отвечающего диссипативным эффектам (эффектам вязкости).

Так как втекание-вытекание газа на границе направлено по нормали, то для оператора гиперболической части на границе систему можно приближенно записать в виде

                                                     ,                                               

где  - вектор консервативных переменных,  - координата по направлению внешней нормали,  - вектор потока.

Систему запишем в виде

.

Здесь  - матрица Якоби, соответствующая потоку , представимая произведением

,

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

Далее систему запишем в виде

                                             ,                                                              

                                               .                                                            

Тогда на границе  граничное условие будет задано следующим образом:

                                          .                                                      

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

4.3                        Расчет нормы вектора вихря

Норма вектора вихря вычисляется по формуле (см. [10]). Для аппроксимации частных производных использован алгоритм, аналогичный описанному в [1]. В каждый сеточный момент времени и в каждой ячейке вычисляется норма вектора вихря. Целью расчета является нахождение этой величины для различных вихревых камер.

4.4                        Пространственная сетка

Ниже на рис. 3, 4 представлен вариант пространственной сетки. Область, в которой будет разыскиваться решение, разбивается на тетраэдры. В той части области, где происходит вдув газа, производится сгущение сетки. На боковых гранях секторной области (рис. 2) строится одинаковое число тетраэдров поворотом относительно оси .

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

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

Из рис. 5 видно, что зависимость нормы вектора вихря от числа Маха близка к линейной. На практике интересны скорости вдува при числах Маха от 0.4 до 0.8. Поэтому в дальнейшем ограничимся расчетом для одного числа Маха  и разных геометрических параметров, предполагая линейную зависимость от  и для них.

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

Ниже приведены  линии тока, построенные в окрестности образования вихря  для сечения камеры плоскостью  (сечение по центру камеры и отверстий вдува).  Линии тока построены для различных моментов времени (после достижения максимума нормы вектора вихря) и для различных геометрических параметров. Угол наклона сечения и количество линий тока выбирались так, чтобы был нагляднее виден вихрь. Оттенки серого показывают  уровень плотности.

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

Из рис. 7-14 видно, что линии тока в окрестности вихря образуют поверхность тора. На рис. 14 виден дополнительный вихрь сверху около отверстия вдува. Его размер зависит от радиуса отверстия. Чем больше это отверстие при фиксированной скорости вдува, тем больше массовый вдув в центр камеры, который разрывает вихри, и тем больше этот дополнительный вихрь сверху около отверстий вдува.

Ниже показаны результаты проведенных трехмерных расчетов в виде линий уровня давления и плотности для разных геометрических параметров и моментов времени после достижения максимума . Единицы измерения: давление - атмосферы, плотность – килограмм на метр кубический. Относительные радиусы отверстия вдува: , , ,  и . По этим графикам можно судить о движении газа в канале. Момент времени выбран так, чтобы поступающий газа занял хотя бы половину объема камеры.

Из рис. 15 – 30 видно, что во всех случаях вблизи отверстия вдува образуется область с пониженным давлением  и плотностью. Из-за этого происходит образование вихря около отверстия вдува. Размер этой области, вообще говоря, зависит от отношения  и скорости вдува (в нашем случае числа Маха). Также в центре камеры образуется зона повышенной плотности (это хорошо видно на рис. 21-22 и 25-26). Размер этой зоны зависит от относительного радиуса вдува  и скорости вдува.

§ 6. Сравнение с результатами двумерной модели

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

Геометрические параметры таковы: относительные радиусы отверстий вдува  (0.03, 0.04, 0.045, 0.05, 0.065), длина отверстия вдува равна 0.8. Число Маха для всех расчетов одинаковое . Начальные и граничные условия такие же, как в трехмерной задаче.

Сравнение рис. 6 с рис. 31 показывает, что двумерная модель качественно верно описывает трехмерную задачу. Более того, она практически так же описывает точку минимума нормы вектора вихря. Можно заметить, что отсутствие третей компоненты скорости в двумерной модели уменьшает норму вектора вихря относительно трехмерного случая (см. рис.6 и рис. 31).

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

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

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

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

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

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

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

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

Из рис. 35 видно, что, если взять радиус вдува слишком большим, то никакого «нормального» (ламинарного) течения быть не может. В зоне отверстия вдува образуются дополнительные вихри.

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

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

Из рис. 36-39 видно качественное и даже количественное (для рис. 38-39) соответствие распределений поля давления в камере для двумерных и трехмерных расчетов. При дальнейшем расчете картины течений и их основные параметры (плотность, давление, скорость) для двумерной и трехмерной задачи, конечно, будут отличаться. Ниже построены линии уровня давления для момента времени после достижения максимума нормы вектора вихря.

Заключение

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

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

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

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

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

Список литературы

1. А.Е. Бутырев, М.П. Галанин, В.Г. Гнеденко, А.В. Переславцев, Е.Б. Савенков, С.С. Тресвятский Математическое моделирование форсунки канала плазматрона в двумерном приближении. // Препринт ИПМ им. М.В. Келдыша РАН. 2007. № 17. 30 с.

2. В.С. Зарубин, Г.Н. Кувыркин Математические модели термомеханики. М.: Физматлит. 2002. – 168 с.

3. И.И. Смульский Аэродинамика и процессы в вихревых камерах. – Новосибирск: ВО «Наука». 1992. – 301 с.

4. Л.Д. Ландау, Е.М. Лифшиц. Теоретическая физика. Т. 6. Гидродинамика. М.: Наука, Физматлит. 1986. – 736 с.

5. Галанин М.П., Грищенко Е.В., Савенков Е.Б., Токарева С.А. Применение RKDG метода для численного решения задач газовой динамики. // Препринт ИПМ им. М.В. Келдыша РАН. 2006. № 52. 31 с.

6. Куликовский А.Г., Погорелов Н.Б., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит. 2001. – 608 с.

7. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction. Berlin: Springer. 1999. –  624 p.

8. R. Eymard, T. Gallouët, R. Herbin. Finite volume approximation of elliptic problems and convergence of an approximate gradient. //Appl. Num. Math. 2001. 37/1-2.  P.p. 31 – 53.

9. В.С. Энгельшт, и др. Математическое моделирование электрической дуги. Фрунзе: Илим. 1983. – 363 с.

10. Марчук Г.И. Методы вычислительной математики М.: Наука. 1977. -  462 с.

11. А. Пуанкаре. Теория Вихрей. Ижевск: НИЦ «Регулярная и хаотическая динамика». 2000. – 162 с.

12. Ильгамов М. А., Гильманов А. Н. Неотражающие граничные условия на границах расчетной области. М.: Физматлит. 2003. – 240 с.