Математическое моделирование течений многокомпонентных сред в кольцевых соплах

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

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

Движение конденсированной фазы в кольцевом сопле, в отличие от ее движения в традиционно применяемом сопле Лаваля, имеет ряд существенных особенностей, определяемых геометрической конфигурацией проточного тракта сопла [1, 2]. Наличие минимального сечения в виде кольцевой щели и резких изменений площади проточного тракта вызывает возникновение значительных градиентов газодинамических параметров, приводящих к существенной неравновесно-сти течения полидисперсного конденсата в кольцевом сопле. Более сложная геометрическая конфигурация кольцевого сопла приводит к значительному искривлению линий тока, как газа, так и частиц, вызывая сепарацию частиц в направлении, поперечном оси сопла. Сепарация, в свою очередь, приводит к интенсивному осаждению частиц из потока на центральное тело и внешнюю обечайку сопла.

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

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

Расчет и конструирование

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

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

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

Система уравнений нестационарного осесимметричного течения многокомпонентной среды в кольцевом сопле имеет вид:

+ V- ур u = 0;

д t дурu „     - дур\

—— + ^V • у р uu + — = r У NC R ( u si ^- u ) ;

д t                    дx

+ ^V ■ у р vu + - yp ^- p = у У NC R ( ^vsi ^- v ) ;

д t                    ду д u р| H - p |

----- ^р + ^V • у р H W = у У N i { Сa i ( 0 Si ^- 0 ) + C R [ u si ( u si ^- U ) + ^vsi ( ^vsi ^- v ) ] } ;

д t p k -1 ( и 2 + v2

= H .

р k I 2 J д у^        i,

_д( + ^V • y ^Ni u si = ^- y ^Ni У k ij Ф _уЩ ;

^д у р si д t

^{+ V} у р si u si = у

^д у р si u si д t

il

Ni У kijФijрsj - рsi У kijФijNj j=1

^{+ V} у р si u si u si = у

l

^Ni C Ri ^{( u - u}si ) ⁺ N i У k j ^р sj ( ^usj ^{- u}si ) ⁺ j = 1

+р si У kij (1 -Ф ij )( usj - usi ) Nj - usi  р si У kij Ф ijNj - Nj У kij Ф ij р si j=1                                          I    j=i                    j=1

^{д у р} s^^{si + V} • у р si v si u si = у [ N i C Ri ( ^{v - v}si ) ⁺ N i У k ij р sj ( ^vsj ^{- v}si ) ^{+ д} t                                                        j = 1

l                                                                                             l                                          i

+Psi E kij (1 - Фij ) (vsj - vsi ) Nj - vsi  PSi E kij ФijNj - Ni E kijФijPSj j=1                                         I j=1                   j=1

^d yc s ° Si P Si d t

^{+ V} y P Si c s ° Si u si

j

= y [N c (° - ° ) + N V k p (E ■ - E ) +

L i ai \          Si у i ij' Si \ j i I i =1

ll

+P si E kij (1 - Ф ij ) Nj (Ej - Ei ) - cs ° si P si E kij Ф ijNj - NiE kij Ф ij p Sj j=1                                            I    j=1

i = 1, -.., l, где uSi =(usi, vSi) - вектор скорости частиц i-й фракции в осесимметричной системе координат;

u = (u,v) - вектор скорости газа; H - энтальпия газа; V = —+ —; Ej -Ei = cS (°Sj -°Si) + dx dy    J’

₊ J u sj

^USi                                                                             ,         XX

■——; c_S - теплоемкость частиц; ° _Si и ° sj - температуры частиц; k jj - коэффициент коагу-

ляции; O j - коэффициент эффективности соударений; N_i - число частиц i- й фазы; C_Ri и c _{a i}

-

коэффициенты взаимодействия между частицами и газом; l - число фракций.

Представленная система уравнений состоит из двух подсистем: подсистема, описывающая движение газа, состоит из пяти уравнений и обозначена (1), подсистема, описывающая движение частиц, также состоит из пяти уравнений и обозначена (2).

В случае учета дробления частиц за счет вращения к представленной системе уравнений необходимо добавить уравнение для квадрата среднего момента количества движения:

^ NM- ^{+ v} „ yN i M i² u_S i + d t

2 B_iM i

( i

l

v 2 _

2ij

i

yN i = y E kN ( M ‘ 2

+

i=j i = 1, ^, l,

-M2У к-ф.NAN;; i ij ij j i j=1

где Mi - моменты вращения частиц i-й и j-й фаз до столкновения; M’j, M'^j - моменты вращения осколков, получающихся при столкновениях частиц i-й и j-й фаз; Bi - момент силы; Ii - момент инерции частицы i-й фазы.

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

Характеристики подсистем (1) и (2) в пространстве { x , y , t } совпадают с траекториями частиц, вдоль которых выполняется столько уравнений совместности, сколько уравнений содержится в этих системах.

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

Граничные условия задаются следующим образом. Для интегрирования системы уравнений на части границы, через которую осуществляется ввод многофазной смеси, задаются энтальпия, энтропия и распределение скорости газа, на стенках сопла использовано условие непротекания газовой фазы v = y ' _xu , на выходной сверхзвуковой границе граничных условий не требуется. В силу характеристических свойств уравнений (1) и (2) граничные условия необходимо задавать только на части границы, через которую осуществляется вдув многофазной смеси. На данной части границы течение предполагалось равновесным (скорости и температуры частиц приравнивались к соответствующим параметрам газа, а момент вращения частиц полагался равным нулю).

Расчет и конструирование

Для коэффициентов Ф у в расчетах использовались два вида функциональных зависимостей:

– экспериментальная зависимость, полученная А.А. Шрайбером и А.М. Подвысоцким, приведенная в работе [4]

Ф у = 1 - 0,247 Re O, ⁴³⁴ Lp-^³³ y -^0,273 ,                                                             (4)

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

35 <  Re j <  385, 5 <  Lp i <  600, 2 < у j <  12;

– обобщенная зависимость, полученная на основе экспериментальных данных и предложенная авторами работы [4]

Ф j = 1,4 + 1,979 X + 0,507 ( 2 X 2 - 1 ) ,                                                          (5)

где X =

( Re

0,572

v ^Re0 J

Lpi

\ -0,153

V ^L P 0 J

A-0,597

^Y ji

Re₀ = 383,6, Lp ₀ = 370,37, y ₀ = 2,37.

Коэффициенты k_ij определялись в соответствии с подходом, принятым в [4].

Для оценки момента распада частиц конденсированной фазы использовался критерий Вебера

(число Вебера) [5]

We = ^{2 P bV} r , ^о в

где р _в - плотность вещества частицы, V - модуль скорости частицы, r - радиус частицы, о _в -коэффициент поверхностного натяжения вещества частиц.

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

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

f ( r ) = —exp V 2 п r ln о

( ln r - In Г 0 ) 2 2 ( ln о ) ²

где r0- математическое ожидание радиуса частиц, r - текущий радиус частицы, ln r0, (ln о)2- математическое ожидание и среднеквадратичное отклонение логарифма радиусов частиц.

Для функции распределения частиц по массам в настоящем исследовании принято соотношение, приведенное в [4]

N ₀ ( m ) = — • exp m

( In m - In m

18 ( In о ) ²

где m – масса частицы, ln σ – среднеквадратичное отклонение логарифма массы, с – константа.

Константа с вычислялась из условия нормировки

l

Ip si i=1

z P

1 - z

где z – массовая доля конден-

сированной фазы в многофазной смеси.

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

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

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

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

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

Соответствующая система уравнений для нестационарного осесимметричного течения в кольцевом сопле имеет вид

l

V • yN ^u = ^- yN i E k ij Ф _у^ ;

j = i

^V" У P si u si = У

il

N i E k ij ^Ф ij P sj ^- P si E k ij ^Ф ij N j j = 1                    j = i

^V У P si u si u si = У

l

^Ni C Ri ( ^{u - u}si ) + N i E k ij ^P si ( ^usj ^{- u}si ) + j = i

l                                                                                               l                                          i

+Psi E kij (1 - Фу ) (usj - usi) Nj - usi PsiE kijФijNj - NiE kijФу Psj j=*                                          V j=i                   j=1

^V" У P si^vsi^usi = У

i

N i C Ri ( v ^- v si ) ⁺ N i E k ij P sj ( v sj ^- v si ) ⁺ j = 1

l                                                                                             l                                          i

+P si E kij (1 - Ф j ) ( vsj - vsi) Nj - vsi P si E kij Ф ij Nj - Ni E kij Ф ij P sj j'=i                                             V    j=i                     j =1

^V У P si c ⁹ si^usi = У

l

Nc^ (9 - 9„-) + NYk,-P„- (E,- -E) + i ai \ si / i ij» si \ j i i j=i

l                                                                           l                                i

+P si E kj (1-Ф ij) Nj( Ej - Ei)-c 9 si P siE kij Ф ijNj — Ni E kj ф i P sj j'=i                                              V i=j                     j=1

^д^ ^У + V- У P u = 0 ; d t

d у p u

d y p v d t

l

^{+ V}- У P ^u l/ ⁺ E p si^usi^usi

i = 1

l

^{+ V}- У P ^{vu +} E p si^Vsi^usi

i = 1

+ ^ P = 0; d x

d УР     _n

^{+ -} p = 0; d x

;

d y p I H - p

V p d t

+ ^V'P Pt< H ⁺ E p si^usi i = 1

2 . 2 A n u_si + V_si c ⁹ _si + si 2 si

= 0;

p k -1Г u2 + v2 — =-- H-- p k V 2 J

Расчет и конструирование

В случае учета вращения частиц к системе уравнений (9) необходимо добавить стационарное уравнение для среднего значения квадрата момента количества движения

V. yNM г^- si

2 B_iM_i ²

⁺ I i

i

yN i = У Z k ij N j ( ^M [ j = 1

; ²

1 ij

+

li

+У ь (1 -ф„ Um'2 -mNi -m2У k,ф,nl

^Ni ^.

ij ij 2ij i j i ij ij j i=j                                                j=1

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

В области сверхзвукового течения стационарные уравнения газовой динамики, полученные из (1), являются уравнениями гиперболического типа. Стационарные уравнения для частиц также имеют только действительные характеристики – траектории частиц. Это позволяет использовать для расчета сверхзвуковых многофазных течений разностные методы, разработанные для решения задач гиперболического типа. Для расчета параметров течения в сверхзвуковой области удобно использовать схему псевдоустановления. В этом случае замена расчетной области на близкую к прямоугольной осуществляется с помощью следующих конформных отображений: L = yj Уо (x), где У о — координаты внешней обечайки, для области движения газа, L i = y^yi (x), где y = у, (x) - уравнение граничной линии тока частиц i-й фракции, для области движения час- тиц. В отличие от дозвуковой части сопла геометрия сверхзвуковой части сопла имеет более «гладкий» характер, поэтому для отображения области течения на расчетную область, близкую к прямоугольной, достаточно рассмотреть, как правило, один вид конформного отображения. После указанной замены переменных система уравнений принимает вид ^—

V-L y o W = 0;

^V- l о Р ^U W '^l v о P -^ ² y 0 y 0 p = ^ y 0 Z N i C Ri ( ^usi ^{- U} ) ; 5 x         5^                  i = 1

- d

^{V -} ^ y о Р vW ⁺ —^ о p - y о p =^ y 0 Z N i C Ri ( v si ^{- v} ) ;

^,L                      i = 1

V - Ly0РHW- = Ly2 Z Ni {Cai (9si - 9) + CRi [Usi (Usi - и) + Vsi (Vsi - v)}] ; i =1

l

V i -L iNi NW i =-L i- y ² N i Z k ij Ф ij N j ;

j = 1

——

i

V- -Ly p W = Ly- N/Zk/Фг/Р^/ iiisiiiiiijijsj

-Р si Z k j ^Ф ij N j

^V i ^{- L} i?i Р si ^ si W i = ^L i y 2 N i C Ri ( ^{и - U}.

l

⁺Р si Z k ij ( 1 ^-Ф ij )( ^usj j = i

-

;

i

^:si ) ⁺ N i Z k ij Р sj ( ^usj j = 1

Usi ) Nj - Usi Р si Z kij Ф ijNj -

- usi) +

i

N j Z k ij ^Ф ij ^Р sj ;

— ⁱ

Vi - LiNiРsivsiW = Lik2 [NiCRi (v - vsi ) + Ni Z kijРsj (vsj - vsi ) + j=1

l                                                  <  l                          i

+Рsi Z kij (1 - Фij ) (vsj - vsi ) Nj - vsi  Рsi Z kijФijNj - Ni Z kijФijРsj j'=i                                             I    j=i                     j=1

v i .4 i y P s c e si W = 4 y

n, Сш (e—e „,) + N ^, k, ,p., (e, — e ) + i a i \ si / i tj'Sjx j ii

J '= 1

l                                                      r l                          i

+P si £ kij (1 — Ф ij) Nj (Ej — Ei) — Cs e si P si £ kij Ф ijNj — Ni £ kij Ф ij P sj

J = ⁱ                                             к     j = ⁱ                   j = 1

^— 2 BM ²

^v i ЛjAAr W + —4 i- y ² = ⁴i y i

Ii

N i £ kN ( M ²

j = 1

^Mi ²

+

li

+У k, (1 — Ф,, M M 2, — M² ) NN, — M 2 N У к, ,Ф, ,- N ,;

ij ij         2ij          i i j i i ij ij j j=ij i = 1, …, l, d —*   d— d —*   d ——      —             — —         —— где v = —k +—j , v = —k +—I , W = uy0k + (v — ^y0u)j , Wi = usiy0ik + (vsi — ^iyiusi)li.

dx    <4        dxd^

Граничные линии тока частиц yi (x) вычисляются из уравнений для траектории частиц каж- дой фракции dyi = vsi dx usi

и начального условия yi (x) = y0i, являющегося результатом численного интегрирования уравнений до-, трансзвукового течения. В случае осаждения частиц на поверхности соплового блока применяется условие «прилипания» частиц к поверхностям соплового блока, а величина yi при- равнивается к значениям y0 , задаваемым уравнением образующей внешней обечайки кольцевого сопла yо =yо(x).

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

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

- = y о 1 (x), - = y 0 2 ( x) , uu где индексами 1 и 2 обозначены соответственно параметры центрального тела и внешней обечайки.

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

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

– геометрические параметры сопла: радиус, на котором расположено минимальное сечение сопла (служит единицей измерения) R = 1; площадь минимального сечения сопла F* = 0,5; радиус

выходного сечения сопла R a = 2,2: угол наклона минимального сечения сопла е * = 45°;

– параметры рабочего тела: давление торможения на входе в сопло p 0 = 120 кг/см² (12 МПа), массовая доля конденсированной фазы z = 0,3;

– распределение частиц по фракциям на входе в сопло задано по нормально-логарифмическому закону.

Рис. 1. Схема кольцевого сопла