Итерационный метод решения одномерных уравнений газовой динамики Текст научной статьи по специальности «Математика»

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

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

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Масштабный переход в химической технологии. Разработка промышленных аппаратов методом гидродинамического моделирования/ Под ред. A.M. Розена. М.: "Химия", 1980, 319 с.

2. Василевский М.В., Мштьцев А.К., Перков В.В., Танков Н.К., Богданов А.Л. Эффективность сепарации частиц в вихревых камерах//Фундаментальные и прикладные проблемы современной механики (доклады Всероссийской научной конференции). Томск: изд-во Томского университета, 1998, С. 276-27?.

3. Справочник по пыле- и золоулавливанию/ Под. обшей ред. A.A. Русанова. - 2е изд. - М.: Энергоатомиздат, 1983. - 312 с.

4. Смит мл. Экспериментальное изучение вихря в циклонном сепараторе. Техническая механика. Труды американского общества инженеров-механиков, русский перевод. 1962, № 4, С. 229-236.

5. Смит мл. Анализ вихревого потока в циклонном сепараторе. Техническая механика. Труды американского общества инженеров-механиков, русский перевод. 1962, № 4, С. 237-248.

6. Аэродинамика закрученной струи. Под ред. Р.Б. Ахмедова. М.: "Энергия", 1977.

7. Абрамов Ю.И., Балусов В.А.. Силин A.B. Некоторые закономерности износа и прогнозирования пылевой эрозии рабочих колес центробежных компрессорных машин. - Тр. ЦКТИ, 1980. Вып. 181. С. 91-100.

8. Шавкин Н.К. Очистка природного газа на магистральных газопроводах. - М., 1973.

9. Потапов О.П., Кропп Л.Д. Батарейные циклоны. - М.: "Энергия". 1977,152 с.

10. Страус В. Промышленная очистка газов. Пер. с англ. - М.: Химия, 1981,616 с.

УДК 533.6.01135

ИТЕРАЦИОННЫЙ МЕТОД РЕШЕНИЯ ОДНОМЕРНЫХ УРАВНЕНИЙ ГАЗОВОЙ ДИНАМИКИ

В.М. ГАЛКИН

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

Введение

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

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

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

Математическая модель

Рассмотрим одномерные стационарные уравнения для идеального совершенного газа [2]:

сЫрЦ

йх

О,

С,

Я

с1(ри2+Р) 1 <М тт,

——-- +--ри~

скх А (Ьс

<Лх А аЬс

у Р и2

=-7—+ Р = рш\ х с |.у |.

у-1р 2

(1)

(2)

(3)

где р, II, Р, Т, Н, Я, у- плотность, скорость, давление, температура, полная энтальпия, газовая постоянная и показатель адиабаты газа, А - площадь поперечного сечения сопла, х - продольная координата, принадлежащая рассматриваемой области С1 иС,-в общем случае нелинейно зависящие от параметров газа правые части уравнений движения (2) и энергии (3). Полагается, что Я и у постоянны, задана площадь А(х), а также известен вид С и Г. Граничные условия заданы на входе в сопло в виде: Я = Н)} и 5 = 5Э, где 5 = Р/ру - энтропийная функция. В качестве априорной информации используется тот факт, что внутри рассматриваемой области число Мбольше нуля, монотонно возрастает от дозвуковой до сверхзвуковой величины и существует только одна точка х*, в которой М(х*) = 1:

М{х„) = \, х„ е(ха;хь);

<1М

сЬс

>0, М(х)> 0, хе[х(д]; М{ха) <\, >\.

(4)

Вместо уравнения неразрывности (1) воспользуемся его интегралом:

рШ=С0, (5)

где С0 - неизвестная константа. Переходя к переменным Н, 5и М перепишем уравнения (2), (3), (5) в следующем виде:

йх

(6)

М = Сп

(7)

(8)

где

р = |-

уу8(2 + (у -\)М )

Р = 8рг, I/ = М

г 2(у - \)Н > 2 +(у -1)М2 у

(9)

Предположим, что для уравнений (6) - (7) в рассматриваемой области известны правые части. Тогда решая задачу Коши для обыкновенных дифференциальных уравнений (6) - (7) с граничными условиями Щ и 50, заданными на входе в сопло, найдем распределение вдоль сопла И и $. Если известна константа С0, и точка х*, в которой М(х*)=1, то из трансцендентного уравнения (8) выбирая дозвуковое или сверхзвуковое решение, вычислим распределение вдоль сопла числа Маха М. Таким образом, имеется система трех уравнений (6) - (8) с пятью неизвестными Н, 5, М, С() и х*. Для решения этой системы ее необходимо замкнуть путем нахождения С0 и х*.

Рассмотрим задачу нахождения С0 и х*. Для этого обратимся к газодинамической функции плотности потока массы [3]:

у+\ р^-!)

Я

2 + (у - \)М'

М.

В точке х* для функции ц и ее первой и второй производной имеем:

щ

м=

йд с!х

й2д

м=

ск

<0.

(10)

Последние два соотношения в (10) являются необходимым и достаточным условием существования в точке М= 1 максимума функции д. А при выполнении условий (4) этот максимум будет единственным. Обозначим через ТУ выражение:

(

М = А

2(7-1) у+1

Я

7+<

2(уЧ) ( } V-

Тогда (8) можно записать в следующем виде:

или

7+1

Щ = С0

7+1

N

/

2 + (у-1)М

Из (12) для функции И, ее первой и второй производной, с учетом (10) получаются следующие соотношения:

(П)

(12)

(13)

С0 еШ __ С0 с1д а1Ы 2Сй Ыд Х С0с12д

д ск д2 с1х ' ск1

Г

йх

\ ил J

д2 ск2

Ы\ = С

(14)

ск

= 0,

й2Ы

ск"

>0.

(15)

Очевидно, что (15) является необходимым и достаточным условием существования в точке М-1 минимума функции И, который при выполнении условий (4) будет единственным. Легко доказать обратное: если при выполнении условий (4) функция #имеет в точке х* минимум, то в этой точке М(х*)=1. С учетом (14) это дает возможность нахождения Со, которое является минимальным значением функции N. а координата положения минимального значения функции Л"дает значение х*:

С0 = 1Шп(7У), х, = х\

шин /V; '

(16)

Отметим, что необходимое условие существования минимума следует из закона обращения воздействия [4], а если в (6) и (7) положить Н=сош и б^соМ, то вышеприведенные выкладки будуг соответствовать экстремальным соотношениям для уравнения Погонио [5, с. 113].

Иногда априорная информация (4) является слишком строгой в том смысле, что могут существовать участки течения с постоянным числом Маха, например, на начальном ци-линдрическом участке сопла с С]=0 и С2=0. Тогда этот участок в соотношениях (16) не рассматривается.

Численный алгоритм метода

Таким образом, для замкнутой системы уравнений (6), (7), (11), (13), (16) с переменными Н, 8, М, N. Сь х* предлагается следующий численный алгоритм:

1. В рассматриваемой области вводится расчетная сетка м и сеточные функции Д, М{, Иь Си, С2), /'=0,1...к, где к - число точек сетки.

2. Задается начальное приближение, например, в виде Н,=Н{], Д =5о.

3. По формуле (11) вычисляется Л/,.

4. Из (16) находится С0 и х*.

5. Из трансцендентного уравнения (13) вычисляется М-„ причем, если х,<х*, то выбирается дозвуковое решение, в противном случае сверхзвуковое.

6. Из (9) вычисляются р, Р, II, а по ним Си и С2, для правых частей уравнений (6) и (7).

7. Вычисляются правые части уравнений (6) и (7).

8. Для обыкновенных дифференциальных уравнений (6) - (7) с граничными условиями Нй, 50 решается задача Коши, которая дает новое приближение Е„ ¿1.

9. При необходимости следующей итерации производится переход на пункт 3.

Для ускорения сходимости в правых частях уравнений (6) - (7) используется нижняя релаксация в следующем виде:

Ь> (17)

где) - номер текущей итерации, ЬИ правая часть, которая использовалась на предыдущей итерации, Ь> - правая часть, которая будет использоваться на текущей итерации, Ъ - правая часть, которая вычислена на текущей итерации, со - параметр, 0<со<1.

Кроме того, более точное значение х* из (16) находится следующим образом. Пусть Ni минимальное значение сеточной функции из (16), тогда применяя из (15) необходимое условие минимума

^ =0

йх г

к параболической интерполяции функции N по трем точкам, получим:

х (х,2-х,+12)Ми +(х1+12-х„12)М1 -Кхи2-х,2)К;+1

А* —

2((хгх;+1)Мы +(х1+гхи)М, Далее используя х* и параболическую интерполяцию можно уточнить С0:

(Х*-Х;)(Х„-Х1 + 1) (Х„-ХЫ)(Х*-Х; + 1) 0»-Х;.,)(Х,-Х;)

Ч -^1-1 т---Г + -^-Г ¡+1-

(хн-х,)(хм-х1+1) ' (хгх;.,)(хгхы) ' (Х1+ГХИ)(Х1+ГХ1)

Тестирование метода

Для описанного алгоритма был проведен ряд расчетов в осесиммет-ричном радиусно коническом сопле, показанном на рис. 1, со следующими параметрами: хиод=4.5, %ход=3.125, Я1~), 91=45°, й=0.625, 02=15°, где все линейные размеры отнесены к радиусу минимального сечения, равного 0.05 м. Число точек сетки £=40. Для решения обыкновенных дифференциальных уравнений (6)-(7) использовалась схема Эйлера второго порядка точности. Уравнение (13) решалось методом Ньютона. Рассматривалось двухфазное монодисперсное течение без учета коагуляции, дробления, испарения и конденсации. Параметры течения были следующими: давление торможения 50-105 Па, температура торможения Г0=3000 К, динамическая вязкость газа р.0=5-105 Пах при 7о, весовая доля второй фазы г=0.4, число Прандтля Рг=0Л, теплоемкость вещества второй фазы с/=1420 ДжДкгК), молекулярный вес смеси 30 кг/кмоль, показатель адиабаты газа у=1.1, плотность вещества частиц рй= 6044 кг/м3, диаметр второй фазы 4=10 5 м. В качестве правых частей уравнений (2)-(3) использовались соотношения из [7]:

с^рми.-и),

Рис. 1. Тестовое сопло. хия)=4.5; утд=3.125; (¡¡=1; 01=45°; [¡2 =0.625; 02=75°

С2= рАсрс,ХГ>'Т)+иЛ(их -Щ),

а в качестве коэффициентов взаимодействия использовались выражения, предложенные Ковано-Дрейком и К. Хен-дерсеном [7, 8]:

6 /л

¿,-р,о Рг

СН — СН\ +С1!2 + С'«3

-Ми: Мы

Мип

Мип - 2 + 0.459 Яе0'55 Рг°33

\ + ЪА2Ми(,

Ее Рг

сВ1 =24^е+5.

... , „_Rej3.65-l.53 Т/7Л 4.33+ ехр(-0.247—)

^ V

1 + 0.353 Т/Т

4.5 + 0.38(0.03 Ке+0.48лУк^) Л1 , А п д , 8

- + 0.1 Л/ +0.2М,

1 + 0.0311е+ 0.48л/кё сю = 0.65, (1 - ехр(- Му /Ке))

\U-Ul \и -II I р 0с1

ехр(-0.5М,/>/кё)

Л 4гЯТ

м

г ^ л071 т

; .V - Ц ^(>.57

где индекс 5 относится к параметрам частиц, ср - теплоемкость газа при постоянном давлении. Параметры частиц находились из стационарных уравнений [6]:

1-

и.

ах

с!Т с и,—= (Т-7Д ёх с$

с граничными условиями на входе в сопло в виде равновесного течения;

их = и; 7у = Т.

Сравнение проводилось с решением, полученным методом установления по явной схеме Маккормака [9] для нестационарного аналога уравнений (1) - (3):

дрА дАрЪт

- + -

дх

■ = 0,

дрЦ д(рЦ +Р) J_(М П2_г ¿я дх А ах

д(рН-Р) | д(рин) | 1 <м с

д1 дх А с1х

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

На рис. 2 показано распределение вдоль сопла числа Маха, сплошная линия результа-ты, полученные предлагаемым итерационным методом, кружочки - методом установления. Для получения гарантированной сходимости в методе установления делалось 500 шагов по времени. В итерационном методе при со=0.3 требовалось 20 итераций. Как видно из графика, результаты графически совпадают.

На рис. 3 приведено распределение коэффициента расхода С газа вдоль сопла:

р11А

С

рЛА,

где индексом * отмечены параметры в сечении Л/(х*)=1, Атт - площадь в минимальном сечении. Сплошная линия соответствует итерационному методу, штриховая линия - методу установления. Так как итерационный метод использует интеграл уравнения неразрывности (13), то расход вдоль сопла получается постоянным, в то время как метод установления дает близкие результаты только около значения М= 1.

Рис. 2. Распределение числа Маха вдоль сопла: 1 - итерационный метод; 2 - метод установления (схема Маккормака}

0.83

-4.5

-3.5 -2.5

-1.5

-0.5

0.5

Рис. 3. Распределение коэффициента расхода газа С вдоль сопла: 7 - итерационный метод; 2 - метод установления (схема Маккормака)

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

И, наконец, влияние параметра о из (17) на скорость сходимости показано на рис. 5. Необходимо отметить, что при малых правых частях уравнений (2) - (3) можно использовать со близкое к единице и небольшое число итераций, большие правые части требуют малое со и большое число итераций. Увеличение со ускоряет сходимость, но увеличивает вероятность появления осцилляции. Поэтому в данных расчетах использовалось со=0.3, что позволяло ограничиваться 20 итерациями, что на порядок быстрей метода установления.

0.850-

0.846

10 20 30 40

50 60 70 80 90 100 Число точек се тки

Рис. 4. Зависимость коэффициента расхода газа б от числа точек сетки. Итерационный метод

20 30

Число итераций

Рис. 5. Влиянием на сходимость: 6- коэффициента расхода газа. Итерационный метод

Заключение

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

СПИСОК ЛИТЕРАТУРЫ

Пирумов У.Г. // Изв. АН СССР. МЖГ - 1967. - № 5. - С. 10-22. Черный Г.Г. Газовая динамика. - М., Наука, 1988. - 424 с. Абрамович Г.И. Прикладная газовая динамика. - М., ГИТТЛ, 1953. - 736 с. Вулис Л .А. Термодинамика газовых потоков. - М, Госэнергоиздат, 1950, - 304 с. Лойцянский Л.Г. Механика жидкости и газа. - М., Наука, 1987. - 840 с.

6. Стерши Л. Е. Основы газодинамики двухфазных течений в соплах. - М., Машиностроение, 1974.

7. Глазунов А.А., Рычков А.Д. // Изв. АН СССР. МЖГ - 1977. -№6.-С. 86-91.

8. Henderson С.В. // AIAA Journal - 1976. - V.14. - № 6. - Р. 707-708.

9. MacCormack R.W. // AIAA Paper. - 1969. - № 69. - Р 354-356.

The iteration method for calculation of stationary one-dimensional flow in the Laval nozzle is proposed. Gas in the nozzle accelerates from subsonic to supersonic speed. The calculations for two-phase flow have shown good agreement with the result, obtained by the pseudotransient method. The iterations number of the of-fered method less by order the pseudotransient method

УДК (519.9+518.5):532

МОДЕЛИРОВАНИЕ ТЕЧЕНИЯ МНОГОКОМПОНЕНТНОЙ УГЛЕВОДОРОДНОЙ ЖИДКОСТИ В ТРУБОПРОВОДАХ

В.В. ЖОЛОБОВ, Е.И. ТАРНОВСКИЙ

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

Для теоретического исследования течения в трубах, как правило, используются модели механики сплошных сред. Физические законы сохранения, примененные к элементарному объему движущегося флюида, приводят к определяющей системе дифференциальных уравнений, гибкость которой обеспечивается введением ряда параметров, идентифицируемых в последствии по эмпирической информации. Неустановившиеся напорные и безнапорные течения в гидравлическом приближении описываются системами дифференциальных уравнений в частных производных гиперболического типа [1]. Уравнения в частных производных этого типа достаточно хорошо изучены и разработан ряд эффективных численных методов для их решения. Привлекательной стороной такого подхода является универсальность и возможность адаптации к реальным условиям. В случае существенного изменения давления при стационарном течении углеводородных смесей и при неустановившихся процессах могут возникать фазовые переходы, приводящие к образованию паровых объемов различного масштаба (от изолированных пузырьков в объеме жидкости до протяженных участков с неполным заполнением поперечного сечения трубы). Аналогичная ситуация возникает при движении паро-газожидкостных смесей. Полости могут исчезать (схлопываться) в том случае, если они заполнены только паровыми компонентами углеводородной жидкости. Различие в механизмах выделения и растворения газовых компонентов приводит к тому, что для растворения газовых компонентов необходимо значительно больше времени. При этом реализуются режимы течения с неполным заполнением трубы при давлениях существенно превышающих давление насыщенных паров жидкости. Для прогнозирования процессов, происходящих при траспортировке углеводородов в этом случае уже необходимо знать неравновесный компонентный состав фаз. Первичный компонентный состав газопаровой фазы зависит от компонентного состава жидкости и может быть определен путем рассмотрения фазового равновесия нефти как многокомпонентной углеводородной смеси. Последующее трансформирование компонентного состава газопаровой фазы необходимо рассматривать в рамках неравновесной схемы фазовых переходов, учитывающей предысторию процесса. Одним из важных этапов построения математической модели, описывающей движение многокомпонентной смеси, является конкретизация функций, задающих интенсивность межфазного переноса компонентов.

Применение принципов термодинамики необратимых процессов [2] позволяет получать феноменологические линейные соотношения для кинетики фазовых переходов [31. При этом коэффициенты, входящие в эти соотношения, зависят от свойств фаз, структуры смеси, геометрии области движения и определяются из эксперимента или детального анализа микродвижений в окрестности поверхностей фазового перехода. Молекулярно-кинетический анализ процессов уноса и осаждения молекул пара на единицу площади межфазной границы приводит к формуле Герца-Кнудсена-Ленгмюра [4,5,6] для предельной скорости фазовых превращений. Использование этого уравнения в инженерных целях требует знания коэффициента аккомодации для составляющих углеводородную смесь компонентов. Относительно этого коэффициента не только не существует теории пригодной для предсказания его величины, но и отсутствуют простые способы ее экспериментального определения. Имеющиеся в настоящее время данные о коэффициенте аккомодации для различных веществ очень противоречивы. В связи с этим имеет смысл рассмотрение упрощенной квазиравновесной схемы на межфазной границе, предполагающей, что несмотря на неравновесность в объеме фаз, на самой границе выполняются условия термодинамического равновесия. Дальнейшим упрощением является принятие более сильного допущения в виде квазиравновесной схемы в объеме непрерывных фаз, предполагающей