Содержание

Введение                                                                                                                             3

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

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

§3.       Выбор функций типа источников – стоков                                                         9

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

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

§6.       Пространственная сетка                                                                                         17

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

Заключение                                                                                                                        28

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

Введение

Существует ряд актуальных проблем охраны окружающей среды (относящихся к экологии): переработка отходов (бытовых, медицинских, радиоактивных), эффективное сжигание низкокалорийных углей, некоторые плазмохимические процессы, а именно, разложение сероводорода, газификация углей и тяжёлых углеводородов. Перечисленные проблемы требуют нагрева газов (в частности, воздуха) до температур около 3000 – 5000ºC для переработки отходов и вплоть до 10000 – 12000ºC для плазмохимических процессов. Эти проблемы могут быть решены с помощью плазменных технологий, использующих плазматроны постоянного тока в качестве нагревателей газа, что вызывает необходимость определения параметров плазматронов на стадии предпроектной проработки.

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

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

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

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

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

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

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

1)    зависимость решения от координаты вдоль оси плазматрона не учитывается; скорость газа вдоль этой же оси равна нулю;

2)    трехмерная задача заменяется двумерной путем добавления к правой части определяющих уравнений газовой динамики функций типа источников - стоков массы, импульса и энергии;

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

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

5)    число  одинаково для всех рассматриваемых форсунок;

6)    при включении стоков полная масса, энергия и момент импульса в сечении канала должны сохраняться;

7)    толщина стенок канала пренебрежимо мала;

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

9)    задача будет решаться в каналах с форсунками, содержащими шесть или восемь отверстий;

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

11) на границе области задана постоянная температура, различная для границы вдува и границы, соответствующей твердым стенкам;

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

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

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

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

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

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

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

1.4           Критерий выбора оптимальной форсунки

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

                                          ,                                       

                                           ,                                         

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

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

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

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

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

Для полной постановки задачи необходимо задать зависимость функций  типа мощностей источников - стоков.

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

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

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

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

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

                                        .                             

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

, , .

Обозначения границ соответствуют рис. 5:  - боковые границы сектора,  - граница, соответствующая стенкам форсунки,  - граница вдува.

На границе  (см. рис. 5) поставлены следующие граничные условия:

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

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

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

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

§3. Выбор функций типа источников – стоков

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

3.1 Уравнение неразрывности (баланс массы)

Правую часть в уравнении неразрывности будем выбирать в виде:

                                                       .                                           

Постоянная  находится из условия сохранения полной массы газа при стационарном движении. Интегрируя уравнение в стационарном случае по области , с учетом граничных условий получим выражение для :

                                                   .                                                

3.2 Уравнение импульса (баланс момента импульса)

Запишем правую часть в уравнении импульса в виде (считаем что ):

                                         ,                                       

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

Запишем закон эволюции момента импульса в цилиндрической системе координат в проекции на ось :

                  .                

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

                                                     .

Интегрируя по области  с учетом достижения стационарности, находим значение константы :

                                            .                                         

Предполагая, что , записываем компоненты функции стока:

где  - угловая координата рассматриваемой точки.

3.3 Уравнение энергии (баланс энергии)

Аналогично стоку массы правую часть в уравнении баланса энергии будем выбирать в виде:

Постоянная  находится из условия сохранения полной энергии при стационарном движении. Интегрируя уравнение по области  с учетом стационарности, получим выражение для :

                                          .                                        

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

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

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

                                  .                                         

                                ,

                                  ,

                                  ,

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

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

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

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

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

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

Правые части: , , , ,

Где   ,

.

Далее обозначение волны «» у безразмерных величин опустим.

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

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

Для численной реализации на базе программного комплекса plasma.m [4] разработан программный комплекс plasma_f.m, адаптированный под данную физическую задачу.

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

5.1.1    Расщепление по физическим процессам

При построении аппроксимации вязких слагаемых и выделения правой части используется прием «расщепления по физическим процессам» (см. [9]).  Запишем задачу в операторном виде:

,

- вектор консервативных переменных, - оператор задачи, - правая часть.

Оператор  представим в виде:

                                                  ,

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

Интегрирование исходной задачи на отрезке  заменяем интегрированием трех задач, которое состоит из следующих этапов:

                                     ,                                   

                                       ,                                   

                                       ,                                  

                                       ,                                                      

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

На этапе решается система уравнений идеальной газовой динамики для декартовой системы координат; на этапе осуществляется учет правой части; на этапе – учет эффектов вязкости.

Все три задачи считаются в отдельно написанных модулях. О методах решения подробно написано в [5, 6]. Расчет потоков в ведется методом HLLC. Аппроксимации правой части описана ниже. Для решения задачи использован алгоритм, основанный на аппроксимации градиента для конечно-объемных схем [7].

5.1.2     Алгоритм аппроксимации слагаемых, соответствующих диссипативным эффектам

Решение подобно решению задачи вида:

Ниже представлен используемый алгоритм её решения (см. [7]).

                                                  ,                                      

                                                     ,

                         .             

Здесь  - шаг по времени; -  объем ячейки; - элемент множества ребер (), соответствующий ячейке K; - поток через ребро ,  - значение в ячейке K(L);  - мера Лебега ребра (в данном случае длина ребра ); () – множество всех внутренних (внешних) ребер, т.е. ребер, не принадлежащих границе области (ребер, принадлежащих только границе области); - расстояние от центра элемента (центра масс треугольника) до ребра.

Cходимость алгоритма доказана в [7]. Недостатком метода является условная устойчивость. К достоинствам следует отнести свойство консервативности.

5.1.3                       Вычисление множителей в правой части

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

Интегралы по контуру вычислялись при помощи того же метода, которым велся расчет потоков в программном комплексе (методом HLLC), а интегралы по области  (всему сектору) заменялись квадратурной формулой:

                                             ,                                           

где - координаты центра i-ой ячейки,  - «объем» i-ой ячейки.  В программе величины  и  вычислялись из соотношений:

                                                    ,                                                 

                                                    .                                                 

В числителе и суммирование ведется  по ребрам границы .  - потоки массы и энергии через j – ое ребро соответственно,  - длина j-ого ребра,  - объем (равный в двумерном случае площади) i - ой ячейки. В знаменателе интеграл по объему заменяется согласно формуле , суммирование ведется по всем ячейкам области.

Такая замена необходима для выполнения дискретных аналогов законов  сохранения полной массы и энергии в области .

Выражение дает для нахождения параметра  формулу, аналогичную и .

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

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

5.3 Другие параметры

Показатель адиабаты в вычислениях взят равным . Шаг по времени . Начальная температура . Начальное давление = 1 атм, давление при подаче газа  атм. Число Рейнольдса Re = 548167. Скорость звука ,  м/c.

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

                                            ,                                          

 - безразмерная плотность,  - для всех форсунок.

Расчеты проведены для форсунок с разным числом отверстий (6 или 8), различных числах Маха (0.2, 0.4, 0.6, 0.8, 1.0), и разных диаметров отверстий вдува. Для форсунки, содержащей шесть отверстий, брались диаметры 2.0, 2.2, 2.5, 3.1 мм, для восьми отверстий брались диаметры 1.0, 1.75 мм при диаметре форсунки 46.0 мм.

§ 6. Пространственная сетка

Ниже на рис. 7, 8 представлены пространственные сетки, построенные с помощью пакета программ PDE Toolbox системы Мatlab7.0 [2].

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

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

   -     в большинстве случаев при построении сетки генератором сеток Matlab возникает проблема с крайним левым треугольником. Выполнение условия периодичности требует написания дополнительной процедуры разбиения этого треугольника;

  -      так как отношение  мало, то возникает необходимость дополнительного сгущения сетки в окрестности области вдува;

  -      во всех случаях сеток рис. 7–10 выполнено равенство .

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

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

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

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

7.2 Распределения параметров газовой динамики в форсунке канала плазматрона

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

Графики линий уровня и векторных полей выполнены с помощью пакета TecPlot 360.0. Среднее число линий уровня выбрано равным 20.  Различное число линей уровня соответствует оптимизации TecPlot для наглядного представления изменения параметров в исходной  области. Все величины на графиках имеют размерность системы СИ (давление – в атм). При рисовании векторных полей показан вектор в каждом узле сетки. Этим вызвано большое число векторов в области вдува (особенно для N = 8), так как там проводилось дополнительное сгущение сетки. Чем большее значение принимает модуль физической величины (в нашем случае - модуль скорости), тем больше длина этого вектора. Все графики построены для момента времени, соответствующего достижению модулем вектора вихря наибольшего значения.

7.3 Обсуждение результатов

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

Из рис. 11 видно, что для форсунок, содержащих шесть отверстий, для чисел Маха 0.2 - 1.0 и рабочего давления на входе атм оптимальным является диаметр отверстий 2.2 мм. При числах Маха 0.8 - 1.0 возможен выбор диаметра отверстий 2.5 мм.

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

Заключение

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

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

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

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

2. Website: http://www.matlab.com.

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

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

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

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

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

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

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

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

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