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

УДК532.517.4: 532.517.2: 532.5.011.1: 544.431.7

К. А. Терещенко, Н. В. Улитин, Р. Р. Набиев

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

СИНТЕЗА БУТИЛКАУЧУКА С ПРИМЕНЕНИЕМ МЕТОДОВ

ВЫЧИСЛИТЕЛЬНОЙ ГИДРОДИНАМИКИ

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

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

Keywords: butyl rubber, fast liquid-phase processes, computational hydrodynamics, Navier - Stokes equation system, copolymerization.

A computational algorithm of hydrodynamics and thermal conditions in the butyl rubber synthesis reactor which is based on standard methods of computational hydrodynamics was offered: equation digitalization was carried out; a combined discrete scheme concept of the obtained system was made by the method of setting on time coordinate and geometric coordinates sweep method.

Введение

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

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

Алгоритм расчета и обсуждение результатов

Исходная система уравнений имеет следующий вид [8-10]. Уравнение неразрывности всей реакционной массы:

—V = 0.

аС'

Уравнения неразрывности по каждой компоненте реакционной массы:

dpY d (pvjY )

f

dt

dC i

pDi +

Mt ScT

л

S2Y

a;

2L+Si-

Уравнения сохранения импульса для сплошной среды в проекциях на оси координат: Я\/к

Р^ + =РРк + ^рк', к = 1,2,3.

Уравнение сохранения энергии для сплошной среды

+ Р\ЧИ = ^ +ткЧ\к + АН, 51 а1

(

Tii = M

avi avj

2 „ Sv, --5iiu—■

3 u a;

, ST

, qi = —,

m SCi

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

^(pK )+-^(pVjK ) = A atv ' a; J ' a;J

+Pk -PopKro-

m +

Mt

JK3

sk

Mt

Jffl3

a; i

Sro

a; i

51 1 ' ас лр ' > ас

Ю _ 2

+азкРк -РзРЮ ,

где р - плотность жидкости; С' - ьтая координата в пространстве; V' - ьтая компонента скорости; Р -

плотность массовой силы; р' - компоненты вектора нормальных напряжений; V' - оператор Гамильтона; рк' = -рдк' + тк' - компоненты тензора напряжений; р - тепловой поток; р - давление; И = СрТ -удельная энтальпия; - символ Кронекера; ц -вязкость смеси; X - теплопроводность смеси;

"т = рК, - турбулентный коэффициент вязкости;

ю

Р0 = 0.09 , а = 5/9 , р0 = 0.075, стК1 = 2, стю1 = 2 - стандартные параметры модели турбулентности

, ЗУ: ..

[8]; Рк =-рУ!'у]'

ЗС;

генерация кинетической

энергии турбулентности; Vj' - пульсация компоненты вектора скорости; верхней чертой обозначены параметры, осредненные по Рейнольдсу; У| - массовые доли компонентов смеси; БОт - турбулентное число Шмидта [8]; Б| - скорость образования ь го компонента; й| - коэффициент диффузии 1-го

компонента; ДИ - выделение тепла вследствие протекания химической реакции.

В общем виде каждое из уравнений системы можно представить в следующем виде

Зи ЗР 50 ЗИ . — + — + — + — = 0

(1)

З Зх Зу Зг где и - вектор параметров процесса (поле импульсов, поле температур, поле концентраций), Р, 0, И - вектора движущих сил, которые действуют на грани элементарного объема, содержащего эти параметры.

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

р " ри " рУ

и = ри , Р = ри2 + Р -^хх , 0 = риУ -Хху

рУи -ТуХ ру2 + р -Туу

где тензор напряжений равен

ти =

"121

Г Зи Зу "1 Зу Зх

Зи Зу Зу Зх

"2 ЗУ

Здесь и в дальнейшем были введены новые упрощенные обозначения по сравнению с исходной системой = и, у2 = V, = х, = у.

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

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

дому параметру, не давая им «исчезать» в результате погрешностей на каждом шаге расчета

ису + +з0 1 СУ = 0.

З^ ^ ^ Зх Зу )

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

Н иСхСу + | (Р

1+1/2,] ' 1-1/2,]

)Су +

+1 (0|,]+1/2 - 0|,]-1/2)Сх = 0, или при переходе к конечным разностям

Ди1^+ + (Р1+1/2,] - Р-1/2,]) + (01,]+1/2 - 01,]-1/2 ) = 0 (2) Д1 Дх Ду

где Д1 - величина шага по времени; Дх , Ду - расстояния между соседними узлами сетки; нижние индексы показывают координаты ячейки сетки, верхние - шаг по времени.

Рис. 1 - Нумерация ячеек в равномерной прямоугольной сетке. Заштрихованная область - малый контрольный объем

Так как поле векторов Б и в в каждой момент времени определяется значением и, то можно связать их приращения во времени с приращениями и (верхний индекс - номер шага по времени):

Рп+1 = Рп +Г—1 5ип+1, 0п+1 = 0п + Г— 1 5ип+1.

^Зи) ^Зи)

Здесь и далее 5 обозначает приращение параметров. Эта схема называется явной схемой решения (движущая сила процесса определяется состоянием параметров на текущем шаге), она не обладает безусловной устойчивостью, то есть решение будет устойчиво только при определенных соотношениях Д1 и Дх. Поэтому лучше использовать неявную (движущая сила на данном шаге зависит от параметров на следующем шаге) или смешенную схему решения [8] (смешанная схема была заложена в алгоритм в дальнейшем):

Р = (1 -а)Рп +аРп+1 = Рп +аГЗр) 5ип+1,

0 = (1 -а)0п + а0п+1 = 0п + ) 5ип+1,

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

ма) решение обладает безусловной устойчивостью [8].

Для удобства моделирования разделим вектор движущих сил на конвективную Рс, Сс и вязкую Ру , Су составляющие

Р = Рс + Ру,

С = Сс + Су,

где

о о

Ру = —ххх ,Су = —хху

-хух —Хуу _

ры " ру

Рс = ры2 + Р ,Сс = рыу

руы ру2 + Р

а также введем более удобный вектор параметров

р " и "

у = ы = ^/Ц

у из/Ц

где Ц , , из - компоненты вектора и.

Тогда вектора Ру и Су можно представить в следующем виде

(3)

Су =

0

аы Зу

-ш — + — ау ах

-ш|2

. ау ау

ар

Зу Зы ау Зу ау

= -ш

ар ах аы ах ау ах

(4)

= М

СХ

ау ах

+м

'СУ

ау ау'

где Мрх, Мру, Мсх , Мсу - матрицы свойств среды.

Возвращаясь к предыдущему вектору параметров, получаем

ру = МРХ

ау аи вя ау аи

+ м.

аи ах ау аи

ру аи ау' ау аи

су = МСХ тл + МСУ тп

аи ах аи ау

где матрица перехода

ау аи

о

и1

о

0

1 и1

ы

Р V

Р

0

1

Р

о

Однако в случае несжимаемой жидкости выражения (3) и (4) могут быть упрощены. При условии

аы ау . — + — = о ах ау

сумма

ару ас

у

ах

будет равна ах,

ау

в проекциях на оси координат

=цАГ 2 аы 1 + 1=

ах ау ах (ах) ау (ау ах)

= ш

= ш

ах,

а

аы 1+ а

ау

ах (ах) ау (ах а

а 2ы

а 2у 1

ах2 ах2

аы ау

—+—

ах (ах ау

а 2ы

а 2у 1

ах2 ах2

= ш

Г а 2ы

а 2у 1

ах2 ах2

ух

ах

уу

ах

а | аы ау 1 а г _ ау , = Ш—| — + — 1 + Ш—| 2— I =

ау ах ( ау ах) ау (ау

= ш

а г ау 1 а г аы

а2ы ау2

а2 у 1

= ш

а2ы

2

а2 у

ау2 = ш

ау (ау) ах (ау у

а Г аы + ау4 + + ау (ах ау) ау2 ау2

То есть без ущерба для точности расчета

( ~2 ~2 \ а2ы а2у

ау2 ау2

вязкие члены можно принять равными

"РХ

о

'аы 1

1

1

ау = МРХ

ах

= -ш

ар ах аы ах ау ах

ау аи

аи ах,

Су =

-ц

= -ц

ф

Зу Зи Зу Зу Зу

ЗУ зи

= — = Мсу——.

ЗУ

— = Моу Зу СУ Зи Зу

Тогда в случае смешанной дискретной схемы выражения для вязких потоков на границах контрольного объёма будут выглядеть так:

+ аМох|^У Т —8ц"*1 =

К' )-,,2 =И)

4-1/2

„п I „п

Р|,] + Р|,]-1 2

„п^п . „п ,,п

Р|,]ии +Р|,]-1ии-1

„п.,п . „п ,,п

Р|,]у|,] +Р|,]-1уи-1

Зи; Зу

-ац

0 0 0 0 1 0 0 0 1

ии-1 + ии

0 2

0 0

2

Р|,]-1 + Р|,]

Р|,]-1 + Р|,] Р|,]-1 + Р|,]

У1 + Уц 0

Р|,]-1 + Р|,]

8Р|,] -8Р|,]-1 Ду

8( ии )-8(Р|,]-1ии-1)

ду

8( )-8(Р|,]-1Уи-1)

ду

Аналогичным образом выражаются остав-

рП+1 сп+1 сп+1

шиеся СУ . .,„ , РУ .,„ ., РУ .,„ . . 4,1+1/2 ' У|+1/2,^ У|—1/2,|

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

Суть схемы заключается в следующем. Если имеется уравнение

Зи Зи „

— + с— = 0, З1 Зх

нужно применять схему 5и| и| - и|-1

- + с-Д1 Дх

8и|. + си|+1 - и|

= 0, при с > 0

Д1 Дх

= 0, при с > 0.

То есть 8и|

■ + с+

и| - и

|-1

+ с _

и|+1- и|

Д1 Дх Дх

с + = тах(с,0), с- = т|п(с,0).

В нашем случае получаем ЗРС ЗС,

= 0,

с = ЗРС Зи + ЗСс Зи

Зх Зу Зи Зх Зи Зу

= ЗГС зу ЗЦ эсс ЗУ Зи

= ЗУ Зи Зх ЗУ Зи Зу Обозначив

ЗРС ЗУ

Аг =

"г

ЗУ Зи

А = ЗСс ЗУ

Сс ЗУ Зи,

в отсутствии вязких потоков, получаем

+ А" -Зи + АС Зи = 0 З1 "с Зх Сс Зу

- + 3"с ЛРс (3"с)

1 Зи

+ 3Сс ЛС

К)

1 Зи

= 0

эи

З1 Гс Гс V Гс' Зх °с °с' Зу где Брс - некоторая вспомогательная матрица (матрица преобразования подобия), позволяющая перейти к диагонализированной матрице Лрс .

Это уравнение мы можем расщепить на 3 отдельных, каждое из которых можно дискретизи-ровать с применением разностей против потока. Применив после этого обратные преобразования, получим

8Ц:

- + 3Р Лг р р

Д1

("с +

+3"с Л"с

+3Сс ЛСГ

+3Сс ЛСс

- К) •К) К)

Дх

1 Ц+ц - и, +

Дх

1 ии - ЦН + Ду

1ЦЧ+1 - и,! = 0 Ду

(5)

где

Гхрс 1+ 0 0 ^ 'яр 0 0 Л

Лр + = рс + 0 ч рс 2+ 0 , Лр - = рс 0 х2с 0

0 V 0 ч рс Л3+у 0 0 V х3с 3 /

Ч = тах(Х,0), Х- = т|п(Х,0),

X - собственные числа матрицы АРс . Аналогичные выражения записываются и для матрицы АСс . Матрица преобразования подобия 3 = ( э2 Б3)

- это матрица, состоящая из компонент собственных векторов Арс (каждый вектор Б| - один столбец).

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

х

и

Р|+1/2,] - Ч ЛРс

Р|-1/2,] : = Ч ЛРс

3|,]+1/2 - З3с ЛС

3|,]-1/2 - З3с ЛС

'и+ Ч ЛР

(Зрс (РС )м/2Л -(РС У-

( Г1 и, + Эзс л3с _ К )Чн + Ззс л3с _(зс )и|,].

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

и— = 1 (3Ц - им), ик = 2 (зц+1 - ц+2),

и - и— + ик

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

и— - Ц, ик - им.

Таким образом, для конвективных потоков получаем

(рс )|+1/2,] -

^-и] зи, + зи^- ип+2,])

- Ар

+Ар

зип - и;-1,]

(( + 3иП + 3и[+и - и]

(зи^- Ц+2,

п ^ (6) +

+аАР

+аАР

( + ип+1,])

2

>и + иП+и)

(би)п+1 +

«г

Аналогичным образом выражаются остав-

шиеся з

п+1 ^п+1

зс

п+1

С|,]+1/2 ' С|+1/2,] ' С|-1/2,] • В конечном итоге, дискретная схема примет

вид

8ип+1 ((рС)|+1/2,] +(РV )|+1/2,] (РС)н1/2,| (РV )н1/2,|)

1-1/2,] V «/|-1/2,]'

А! Ах

((3С )|+1/2, | + (3V )|+1/2, -(3С )И/,] -(3V )И/,])

Ау

- о,

где

(РС )|+1/2,] -(РС )|+1/2,] +

(АР+) (биГ1 +(АР ) (5и)"+1 V рс + / |+1/2,Л А] V Рс-/|+1/2,Л

|-1/2,]

+а

(Арсх^аий+к-ц,«

К Г'] (3С )|,.+1/2 -(3С )"

/|,]+1/2

+а

+а

К + ))+1/2 «Г +(А3с-))+1/2 (аи)-(3С )|,]—1/2 -(3С £-1/2 +

(А3с+) «1+(А3с-и н; (Рv )|+1/2,] - (V )п+1/2,] + аЦ+1/2,] "" [(би)п++11 -(8и); (Рv)|-1/2,. -(р/)м/2,|+аьп-1/2,]^[К+1 -«

п+1

(3v )|,.+1/2 -(3^) +1/2 +аМп]+1/2 -УМ+1- (би)|

(3V )|,]—1/2 -(^ ))-1/2 +аМ.1/2

зи

ЗV ^

13X1 ай|

I - единичная матрица.

АУ

_1_ АУ

|,]+1 \п+1

г

— - Мрх| — I, N - М,

После преобразований выражения (6) полу-

чаем

бип]+1[|^((Ар+)п -(Ар )п |,] 1 АХ ^ Рс+^+1/2,1 V Рс )|-

|+1/2,]

/|-1/2,]

((А3+)п -(А3 )п Ау (V 3с + А,]+1/2 V 3с /и-

|,]—1/2

аА! Ах2

+5и|

(| п ,п ) аА (п Nn )1

у'—1+1/2,] —1-1/2] "Ау2 V |,]+1/2 'ЛН^-Т

п+1 ( аА! |+1,]

(Ар )п

V Рс /|+1

. аА^п

Ах у 'с ;|+1/2,] АХ аА^ \п аА! Ах

ГН+1/2,]

+8ип-11 К п

(

+би|

п+1 |,]+1

аА!

д7

К

Ч-1^] Ах' аА!

Л

|,]+1/2 Ау

N

2 |,]+1/2

+5и|,

См

(А3с+)и-1/^

А!

((РС )+1/2,] (РС )-1/2,] + )!+1/2,] )-1/2,])

Ах

А!

((3С )п]+1/2 -(3С )п]-1/2 +(3V )п]+1/2 -(3У )п]-1/2 ) АУ ,

или в общем виде (все известные на шаге по времени параметры объединены в матрицы А| ], В| ], С| ],

°| ,], Е| ,], Аип]):

с

А„ 5У-1 + Би 5УП++11] + Ои бип.^^1 + Ри 5^]+ +

I,] I,]

.Е^и]1, -ди,.

(7)

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

Ац бУ п + Ви бУ 1.1,] + Ои 5УП_1,] -Д^, (8) А^1 + Ц/и"^ + Б^и^ = Д1й П], (9)

где верхней волнистой линией показаны величины, получаемые на промежуточном этапе шага по времени.

После этого каждая из систем (8) и (9) решается методом прогонки по координатам [8-10]. Система

А|,] 5и Щ + В|,] 5и ^ + О|,] 5и П_1,] -ди"

переписывается к виду

А5и" = Р",

Где

5и" =

" 5и"] " " Ди"] "

5и 2,] ,рп = Ди2,]

5и :-1,] ду:-1,]

_ 5и:,] _ _ ди:,] _

А =

А. В.

2,]

А2,] В3,] .. 0 0 0

0 О2,] А3,] .. 0 0 0

0 0 0 .. Ат-2,] Вт-1,] 0

0 0 0 .. От-2,] Ат-1,] Вт,]

0 0 0 .. 0 От-1,]

Матрица А приводится к наддиагональному виду путем следующих преобразований

1-1,]'

А'1, ] - А1,] ,

О| :

= А| : 1,] В|

1 ,] А'

А' =

А'г] В2,] 0 ... 0 0 0

0 А'2,] В3,] ... 0 0 0

0 0 А'з,] ... 0 0 0

0 0 0 .. А т-2,] Вт-1,] 0

0 0 0 ... 0 А т-1] Вт,]

0 0 0 ... 0 0 А т,]

Соответственно

ди

1]

ди'2,]

ди':-1,] ди':,]

где Ди'1; =Ди1], Ди'! -ди" -

А'

1-1,]

-Ди'!-1,].

После этого 5ип, начиная с 5ит ], можно вычислить по формулам

5и

т,] А'

5и П =■

Д и'П: - В| ;5иП.

I,]

I,]

1+1,]

А'

А ¡,]

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

Заключение

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

Р

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

Литература

1. Улитин Н.В. Кинетика катионной сополимеризации изобутилена с изопреном в неизотермической постановке задачи / Н.В. Улитин, К.А. Терещенко, Р.Р. Набиев, Т.Р. Дебердеев, Р.Я. Дебердеев, Э.Р. Гиззатова, С.И. Спивак // Вестник Казанского технологического университета. - 2013. - Т. 16, №19. - С. 193-200.

2. Ибрагимли Д.Ш. Идентификация процесса полимеризации при разработке АСУ ТП производства бутилкау-чука / Известия вузов. Нефть и газ. - 1985. -№7. - С. 8690.

3. Сангалов Ю.А. Полимеры и сополимеры изобутилена: Фундаментальные проблемы и прикладные аспекты. Уфа: Гилем, 2001. 384 с.

4. Маркина Е.А. Синтез бутилкаучука с использованием модифицированной каталитической системы на основе хлористого алюминия: автореф. дис. канд. хим. наук: 18.03.2010, 21.04.2010 / Маркина. - Казань, 2010. С.19.

5. Yasuyuki Tanaka // Progress in Polymer Science. 1989. V. 14. Iss. 3. P. 339-371.

6. Kennedy J. P. // Journal of Macromolecular Science: Part A:Pure and Applied Chemistry. 1982. Vol. A18. Iss. 1. P. 129-152.

7. Кеннеди Дж. Катионная полимеризация. Пер. с англ. под ред. П. Плеша. Москва: Мир, 1966. 584 с.

8. Аникеев А. А. Основы вычислительного те-плообмена и гидродинамики / А. А. Аникеев, А. М. Молчанов, Д. С. Янышев. - М. : URSS, 2009. - 149 с.

9. Франк А. М. Дискретные модели несжимаемой жидкости / А. М. Франк. - М. : Физматлит, 2001. - 206 с.

10. Труды ИБРАЭ. Вып. 3. Методы вычислительной гидродинамики для анализа безопасности объектов ТЭК / науч. ред. В. Ф. Стрижов. - М. : Наука, 2008. - 207 с.

11. Захаров В. П. Физико-химические основы протекания быстрых жидкофазных процессов / В.П. Захаров, А.А.Берлин, Ю.Б.Монаков, Р.Я. Дебердеев. - М: Наука, 2008. - 348 с.

© К. А. Терещенко - асп. каф. технологии переработки полимеров и композиционных материалов КНИТУ, nucleurmind@yandex.ru; Н. В. Улитин - д.х.н., проф. каф. технологии переработки полимеров и композиционных материалов КНИТУ, n.v.ulitin@mail.ru; Р. Р. Набиев - асп. той же кафедры, nabievrafit@mail.ru.