CC BY
32
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XXV 1994 №1-2
УДК 532.529
МЕТОД РАСЧЕТА МЕДЛЕННЫХ СТАЦИОНАРНЫХ ДВУХФАЗНЫХ ТЕЧЕНИЙ С БОЛЬШИМ ОБЪЕМНЫМ СОДЕРЖАНИЕМ ЧАСТИЦ
А. С. Войновский, М. Е. Зайцев
Разработана псшунеявная конечно-разностная схема, эффективная для расчета медленных течений двухфазной смеси, позволяющая рассчитывать течения с большим содержанием дисперсной фазы и имеющая более слабое ограничение на шаг интегрирования по времени, чем явные схемы. Схема имеет второй порядок по пространственным переменным на установившемся решении и использует метод расщепления по физическим процессам на основе разделения характерных времен.
В настоящее время для расчета двухфазных стационарных течений — смеси газа с твердыми частицами на основе установления по времени — разработано множество численных методов (см. [1—3] и др.). Большинство из них построено на основе использования явных конечно-разностных схем. Однако когда в течении характерные скорости меняются на порядки, использование явных конечно-разностных схем практически невозможно из-за жестких ограничений на шаг интегрирования по времени. В связи с этим в данной статье предлагается полунеявная схема для расчета двухфазных течений, имеющая более слабое ограничение на шаг интегрирования по времени.
1. Рассматривается нестационарное течение двухфазной смеси — газовзвеси. Уравнения, описывающие течение, основаны на физических законах сохранения массы, импульса и энергии, выписанных отдельно для каждой фазы с применением гипотезы о взаимопроникающих континуумах [1], и используют следующие допущения:
— частицы — сферы одинакового диаметра г/, не взаимодействующие между собой,
— температура вещества частиц во всем их объеме постоянна,
— несущая фаза — идеальный газ, вязкость, теплопроводность учитываются только при межфазном обмене импульсом и энергией,
— вращение частиц и межфазный массовый обмен отсутствуют,
— силой Бассе и эффектом, связанным с присоединенной массой, при расчете межфазного силового обмена пренебрегается.
Тогда уравнения, описывающие совместное движение газовой (индекс 1) и дисперсной (индекс 2) сред, можно записать в следующем виде:
%- + с11у(р1К1) = 0, (1)
01
+ вшКрЛ • УО + (1 - а)8Пйр + /а = 0, (2)
ах
+ <Му(Р1ВД = а Г28гас\р - /12У2 -д-
^[/»(а^+а-о)^)], (3)
% + с1Мр2^2) = 0, (4)
01
д(-р^ + grad(p2F2 • У2) + а упй р - /12 = 0, (5)
8^г) + ^(р2е2^) = Я,
Е1=Р(У^1Ж.
Р1(ее - 1) 2
(6)
Здесь а — объемная концентрация дисперсной фазы, р, р, У, е, Т— давление, плотность, вектор скорости, внутренняя энергия и температура. Силу межфазного взаимодействия /12 и тепловой поток д можно записать в следующем виде [1]:
/12 = ^¥-р№ - Уг)\У1 ~ Г2| - Уг\
уе2
4(1 г.
" р( 1 - «)
6а ИиХ „ . „
д = -~Г (71 - Т2) = а В
Р1(® -1)
Здесь у = с„ / с2, ее = ср / си, ср> сц — удельные теплоемкости газа при постоянном давлении и объеме, с2 — удельная, теплоемкость вещества частиц, X —теплопроводность газа.
Выражения для са, N11 для различных режимов течения с учетом эффекта «стесненности» потока при больших значениях а берутся из [1]. Уравнения (1)—(3) описывают движение газовой фазы, (4)—(6) — дисперсной.
2. Для численного решения системы (1)—(6) используется метод расщепления по координатам и физическим процессам на основе разделения характерных времен. Рассмотрим уравнения (1)—(3) при а = 0 для течения чистого газа, которые описывают:
1) конвективный перенос массы, импульса и энергии,
2) воздействие на поток градиента давления и дивергенции скорости.
При численном интегрировании данных уравнений по явной конечно-разностной схеме условия устойчивости для расчета конвективного переноса имеют вид
т <сД-, с = const, (7)
м
а для расчета второго процесса
Т < Сг—.-, с = const. (8)
\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)(i»j - v2)
(1 - a) " (1 - a) ’
^_^Щ-а)£ + аВ1е dt dx pj p2 (ае -1) P2
(9)
in
В плоскости переменных (х, /) вводится в рассмотрение разностная сетка. Изменение величины во времени характеризуется индексом и, а по координате — /. Рассмотрим алгоритм определения значения параметров на (и+1)-м слое по известным величинам на и-м слое. На первом шаге рассчитывается по явной схеме следующая система уравнений:
^. + 0 5^- = 0 а/+и,эах и’
др др Л
др, др.
« + ^1 -г-®- =
а/ ^ дх ’
^. + 0,5 **
01
дрг
дt
+ ■02
дх
М, = 0,
ах ар2 _
ае, эс2 „
а/ дх
(10)
с начальными условиями на и-м слое и определяются промежуточные значения параметров р, р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>
дР\ОВ аВу (ае - 1) ^ _ ар! А(гг - 1) (о, - у2)
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 )_
(14)
Сі=-
1
Со =
Р( 1 ~ «)
аВ(1/рі + у / р2) [_ (ае - 1) р2
(ае-1)_________ І ¿^(¿і - ¿2)2
е2
(1 - а) (1 / р! + у / р2)^
^Рі Рг
р{\ - а)
Рі(ае - 1)
-уе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.
(15)
На рис. 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
