Метод расчета медленных стационарных двухфазных течений с большим объемным содержанием частиц

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

Том XXV 1994 №1-2

УДК 532.529

МЕТОД РАСЧЕТА МЕДЛЕННЫХ СТАЦИОНАРНЫХ ДВУХФАЗНЫХ ТЕЧЕНИЙ С БОЛЬШИМ ОБЪЕМНЫМ СОДЕРЖАНИЕМ ЧАСТИЦ

А. С. Войновский, М. Е. Зайцев

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

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

1. Рассматривается нестационарное течение двухфазной смеси — газовзвеси. Уравнения, описывающие течение, основаны на физических законах сохранения массы, импульса и энергии, выписанных отдельно для каждой фазы с применением гипотезы о взаимопроникающих континуумах [1], и используют следующие допущения:

— частицы — сферы одинакового диаметра г/, не взаимодействующие между собой,

— температура вещества частиц во всем их объеме постоянна,

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

— вращение частиц и межфазный массовый обмен отсутствуют,

— силой Бассе и эффектом, связанным с присоединенной массой, при расчете межфазного силового обмена пренебрегается.

Тогда уравнения, описывающие совместное движение газовой (индекс 1) и дисперсной (индекс 2) сред, можно записать в следующем виде:

%- + с11у

01

+ вшКрЛ • УО +

ах

+ <Му(Р1ВД = а Г28гас\р - /12У2 -д-

^[/»

% + с1Мр2^2) = 0,

01

д

8^г) + ^(р2е2^) = Я,

Е1=Р(У^1Ж.

Р1(ее - 1) 2

Здесь а — объемная концентрация дисперсной фазы, р, р, У, е, Т— давление, плотность, вектор скорости, внутренняя энергия и температура. Силу межфазного взаимодействия /12 и тепловой поток д можно записать в следующем виде [1]:

/12 = ^¥-р№ - Уг)\У1 ~ Г2| - Уг\

уе2

4(1 г.

" р

6а ИиХ „ . „

д = -~Г (71 - Т2) = а В

Р1

Здесь у = с„ / с2, ее = ср / си, ср> сц — удельные теплоемкости газа при постоянном давлении и объеме, с2 — удельная, теплоемкость вещества частиц, X —теплопроводность газа.

Выражения для са, N11 для различных режимов течения с учетом эффекта «стесненности» потока при больших значениях а берутся из [1]. Уравнения (1)—(3) описывают движение газовой фазы, (4)—(6) — дисперсной.

2. Для численного решения системы (1)—(6) используется метод расщепления по координатам и физическим процессам на основе разделения характерных времен. Рассмотрим уравнения (1)—(3) при а = 0 для течения чистого газа, которые описывают:

1) конвективный перенос массы, импульса и энергии,

2) воздействие на поток градиента давления и дивергенции скорости.

При численном интегрировании данных уравнений по явной конечно-разностной схеме условия устойчивости для расчета конвективного переноса имеют вид

т <сД-, с = const,

м

а для расчета второго процесса

Т < Сг—.-, с = const.

\v\ + a

Здесь v — скорость потока, а — скорость звука, т, А — шаги разностной сетки по времени и пространству. Отсюда следует, что для повышения устойчивости всей схемы — увеличения максимально допустимого шага интегрирования, ограниченного условием (8), до значения т, ограниченного неравенством (7), необходимо рассчитывать распространение возмущений из-за градиента давления по неявной абсолютно устойчивой схеме либо по схеме, устойчивой для т из (7). Данный подход используется и при расчете течения двухфазной смеси. Отметим, что при |и|« а условие (7) дает существенно более слабое

ограничение на г, чем (8). Это обусловливает эффективность расчета с

т из (7) для медленных течений. Численное решение уравнений (1)—(6) на каждом элементарном шаге по t расщепляется на три этапа: на первом проводится расчет конвективного переноса по явной конечно-разностной схеме, на втором — по неявной схеме рассчитывается изменение параметров потока за счет градиента давления и дивергенции скоростей частиц и газа, на заключительном шаге рассчитывается межфазное силовое и энергетическое взаимодействие по точным аналитическим формулам.

Преобразуем уравнения (1)—(6) к виду, удобному для интегрирования. Для простоты рассмотрим одномерное течение. Тогда, используя уравнение состояния, формулы для межфазного взаимодействия и считая, что плотность материала частиц р°2 = р2 / а = = const, можно записать законы сохранения импульса и энергии для обеих фаз в виде

^L + 4*1+<L^> £ + аЛ(ч-ч)-.0,

dt 1 dx р, дх 1

dv2 dv2 a dp \ n

—- + th—- +--------— - a A(Zh - Lh ) = О,

dt 1 дх p2 дх

dp dp aер d

— + V\— +-----------------------

dt dx (1 - a) dx

olB

f(l - a)t>i + a021 +--------------------P ~

Pi

aByeziae - 1) _ a.dpjiae - 1)

(1 - a) " (1 - a) ’

^_^Щ-а)£ + аВ1е dt dx pj p2 (ае -1) P2

in

В плоскости переменных (х, /) вводится в рассмотрение разностная сетка. Изменение величины во времени характеризуется индексом и, а по координате — /. Рассмотрим алгоритм определения значения параметров на (и+1)-м слое по известным величинам на и-м слое. На первом шаге рассчитывается по явной схеме следующая система уравнений:

^. + 0 5^- = 0 а/+и,эах и’

др др Л

др, др.

« + ^1 -г-- =

а/ ^ дх ’

^. + 0,5 **

01

дрг

дt

+ 02

дх

М, = 0,

ах ар2 _

ае, эс2 „

а/ дх

с начальными условиями на и-м слое и определяются промежуточные значения параметров р, р1; р2, ¿2, являющиеся начальными

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

дуг | (1 - а) др = 0

а/ р[ & ’

*И + ±.эр=0> а/ р2 ах

~[(1-а)^+а1>2] = 0,

др гер д/ (1 - а) Зх

^й+Р2«а = 0, ^+Р1«а.=о

а/ дх дt 1 ах

или в векторном виде

|£ + П/ = 0, / = («1, 02, /», Р1, Рг)1.

Записывая уравнения (11) в конечно-разностном виде и проводя факторизацию аналогично [4], для решения системы (11) можно использовать неявную схему расщепления [4] с несогласованным стабилизирующим оператором С^1) с аппроксимацией производных с первым порядком точности:

сыЫ=-ф)~/,

X

С(1) = (Е + сот + «ВТ п?>>.

О = + П2 •

Здесь ш — параметр схемы, Е — единичная матрица, операторы Ох и П2 имеют вид:

П2 =

0 0 (1-а) д Р1 дх 0 0

0 0 а д Р2 дх 0 0

«1 = вер д ааер д п дх (1 - а) дх 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

о 0 0 п еер(х>2 -Оі) д

У/ и и V (1-а) дх

д 0 0 0 0

0 д Р2~ дх 0 0 0

оФ— разностные аналоги операторов Пу, аппроксимирующие их с первым порядком точности, П(2) аппроксимирует О со вторым порядком точности с помощью направленных разностей. Схема (12) в дробных шагах имеет вид:

{е + атО^)^^= -тП^/,

Г-}*4т-

ш

В результате расчета на данном шаге определяются значения плотностей газа и частиц на (и+1)-м слое р”+1 = р^, р2+1 = р2, а"+1 = а и промежуточные значения параметров £^, £¡2, р. На завершающем этапе расчета параметров на (и+1)-м слое, учитывающем влияние меж-фазных силового и энергетического обмена, решается следующая система уравнений:

+ аА(г\ - = О,

- а.А(х\ - ъ) = О,

2 (13>

дР\ОВ аВу

5/ рі (1-а) 2 (1 - а) ’

де-> аВ(\ - а) аВу п

——--------- ----— р +--------е2 = О

ді РіРг(ае -1) Р2 .

с начальными условиями ^(0) = ^(О) = ¿2, р(0) = р, е2(0) = ё2> в

результате решения которой определяются параметры на (и+1)-м слое по Решение (13) имеет вид:

у1"+1 = 0,5(і\ + ог) + 0,5(^ - ¿¡2)ехр(-2ат./4),

Ц2И+1 = 0,5^ + 02) - 0,5(5|. " Уг)ехр(-2«ітА), ^1 = а^1)С1+С2ехр

(ае - і)о.^4(с?1 - ¿і;)2 [атДур^І / рі + у / р2) + р2] р25(1-а)(1/рі +у/рг)2 ’

^+, = «3^ _ ¿^аехр[_іті(і /-Л*1/п)\-

аА{щ ~ Уг)2[-тВ{\ / р! + у / рг)]

Р2г(і/Рі+у/р2)2

Константы Сі, С2 определяются следующим образом:

-ат В

1 Рі Р2 )_

Сі=-

1

Со =

Р

аВ(1/рі + у / р2) [_ (ае - 1) р2

(ае-1)_________ І ¿^(¿і - ¿2)2

е2

(1 - а) (1 / р! + у / р2)^

^Рі Рг

р{\ - а)

Рі

-уе2

При решении двумерной задачи необходимо предварительно провести расщепление по координатам на одномерные подзадачи с производными по координатным осям [4].

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

+

2

Эх,.

Производная от потока (с^) в узле / аппроксимируется на шаблоне (/ - 2, / -1, /, / + 1) следующим образом. Сначала анализируется знак (щ)”. Пусть для определенности (ц)” >0, тогда через точки функции

потока (1^)/_2, (в?)/-1, (*£)/ И (|£)/_1, (ф/, (^Хч! проводятся две параболы и определяются их кривизны и к2 соответственно. При к\ • к2 < 0

ди}

дх

В случае совпадения знаков кривизн выбирается тот интервал, на котором одна из этих парабол имеет минимальную по модулю величину кривизны, и подсчитывается односторонняя разность на этом интервале. Например, при |&1|<|&2| производная аппроксимируется следующим образом:

д(°12)/ (р?)<-2 ~ 4(°2)м + 3(у}),

дх 2А

Если (с^),- < 0, то используются узлы сетки (/' - 1, /, / + 1, / + 2). После определения ^ аналогично определяются остальные параметры газовой фазы. Затем таким же образом определяются параметры дисперсной фазы £¡2, Р2, ¿2-

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

При определении стационарного решения уравнений (1)—(6) на основе их численного интегрирования установлением по / порядок аппроксимации определяется порядком аппроксимации операторов правой части (12) и равен 2. Исследование на устойчивость линейного

аналога уравнений (11) методом Фурье показывает, что схема (12) абсолютно устойчива при со £ 1. Отметим, что при использовании согласованного оператора С в левой части (12) (т. е. порядок аппроксимации С равен 2 и порядку аппроксимации П) схема устойчива в линейном приближении при со £ 0,5. Однако в этом случае система уравнений получается пятидиагональной. При реализации граничных условий для приращений £ при решении (12) используется алгоритм из работы [5], обеспечивающий их абсолютную устойчивость. Заметим, что в двумерном случае даже при со > 1 абсолютная устойчивость схем расщепления с несогласованным стабилизирующим оператором может нарушаться [4]. Однако и при этом, как показывают расчеты, условие вида (7) обеспечивает устойчивость расчета на этом шаге. Неявная аппроксимация источниковых членов при расчете межфазного обмена на третьем этапе обеспечивает устойчивый счет при приближении к равновесному течению двухфазной смеси.

4. Методические исследования разработанного алгоритма были проведены при решении задачи о распаде разрыва в двухфазной смеси — воздух и частицы кварцевого песка, теплофизические параметры которых приведены в [1, 6]. Диаметр частиц <1 = 60 мкм. В качестве характерных параметров выбирались параметры среды низкого давления. При / = 0 задается следующий разрыв параметров:

Р — А) = 1 • Ю Па, Тх = Т2 = То = 293 К, щ = 02 =: 0»

• . .. кг _

р, = 1,21 —а = а<) при х < 0, мЛ

р = 5ро, 7! = Т2 = 293 К, ^ = иг = 0» р\ = 6,05

м-*

а = а0 при х > 0.

На рис. 1 приведены результаты сравнения численного решения данной задачи по разработанному алгоритму (линии I) при а0 = 0,001 в момент времени / = 0,01464 с с расчетом из работы [6] (линии 2), проведенным модифицированным методом крупных частиц [7]. Показаны распределения безразмерных температур газа (сплошные линии и кружки) и частиц (штриховые и крестики). Число узлов сетки 7У=361, для полунеявной схемы время счета — 2 мин на РС АТ 386,

К = т^° + ^ = 1,2. Расхождение с результатом из работы [61 составляет И

не более ~ 5 %. Время счета данной задачи в работе [6] составляло 10—15 мин на БЭСМ-6. На рис. 2, а, б приведены распределения давления и относительной объемной доли частиц а / а0, полученные при расчете распада разрыва (15) с а0 = 0,05 в моменты времени / = 0,59 • 10_3 с (линии /) и / = 0,118 10-2 с (линии 2). В этом расчете число Куранта К = 2,0, N = 101. При К = 3, когда нарушается условие (7), устойчивое решение получить не удается.

л.

о

_____________I

х,м 5

Рис. 1. Результаты сравнения численного решения данной задачи по разработанному алгоритму при а0 = 0,001 в момент времени і = 0, Ь1464

Рис. 2. Распределения давления и относительной объемной долр частиц а / а0

В качестве двумерной тестовой задачи для расчета методом установления рассматривается истечение газовзвеси (воздух и частицы углерода) с большим содержанием дисперсной фазы из осесимметричного канала генератора дисперсных частиц в вакуум. Параметры потока на срезе канала генератора выбирались следующими: ра = 0,2 • 105 Па, аа = 0,-6, Т1а = Т2а = 300 К, ии = и^а = 2 м / с, х\а =

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

Рис. 4. Зависимость отношения времен счета на ЭВМ по схеме из работы [8] %2 и измененной схемы X! до заданного времени с начала момента истечения

= ^ = 0 (и, V — проекции V на оси х, г цилиндрической системы координат), радиус канала Ла = 5 см.

Расчеты проводились на основе изложенной схемы (линии 1 на рис. 3) и методом [8] (линии 2), в котором параметры газовой фазы рассчитывались методом Годунова с модификацией Колгана [9], а параметры частиц — с помощью лагранжева этапа метода крупных частиц [7]. На рис. 3 приведены положения сепаратрис дисперсных частиц и линии равных скоростей частиц в установившемся решении. Максимум расхождения результатов расчета по двум схемам наблюдается для скорости газовой фазы и не превосходит 5%. Для установления течения в рассматриваемой области достаточно 0,15 с.

Сравнения скоростей счета по двум схемам иллюстрирует рис. 4, на котором приведена зависимость отношения времен счета на ЭВМ по схеме из работы [8] т2 и изложенной здесь схемы до заданного времени с начала момента истечения. Видно, что сначала, когда в расчетной области преобладают высокие скорости потока, использование изложенной схемы менее выгодно, чем схемы [8]. Затем, когда область течения с малыми скоростями начинает занимать большую часть расчетной области, изложенная схема обеспечивает большую скорость счета. В итоге выигрыш по времени счета для получения установившегося решения по сравнению с явной схемой [8] составляет -60%. ’

Авторы выражают благодарность И. Э. Иванову за помощь в работе.

1. Нигматулин Р. И. Динамика многофазных сред. — М.: Наука, 1987.

2. Крайко A. H., Нигматулин Р. И., Старков В. К.,

Стернин JI. Е. Механика многофазных сред // Итоги науки и

техники. Гидромеханика. Т. 6. — М.: ВИНИТИ, 1972.

3. Гилинский М. М., Стасенко А. Л. Сверхзвуковые газодисперсные струи. — М.: Машиностроение, 1990.

4. Ковеня В. М., Яненко H. Н. Метод расщепления в задачах газовой динамики. — Новосибирск Наука, 1981.

5. Войновский А. С. Об одном варианте метода расщепления и неявной реализации граничных условий для решения уравнений Навье—Стокса в криволинейной системе координат // Ж. вычисл. матем. и матем. физ.—1990. Т. 29, № 9.

6. Губайдуллин А. А., Ивандаев А. И., Нигматулин Р. И. Некоторые результаты численного исследования нестационарных волн в газовзвеси // Изв. АН СССР, МЖГ.—1976, № 5.

7. Белоцерковский О. М., Давыдов Ю. М. Метод крупных частиц в газовой динамике. — М.: Наука, 1982.

8. Войновский А. С., Зайцев М. Е. Нестационарное истечение в вакуум газовзвеси с большим содержанием дисперсной фазы // Изв. РАН, МЖГ.—1992, № 2.

9- Колган В. П. Применение принципа минимальных значений производной к построению конечно-разностных схем для расчета разрывных решений газовой динамики // Ученые записки ЦАГИ.—1972. Т. 3, № 3.

Рукопись поступила 18/V1992