УДК 517.958:531.12:66.021.1
Ч. Х. Исянов, Ф. Г. Ахмадиев, М. И. Фарахов,
И. Г. Бекбулатов
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА РАЗДЕЛЕНИЯ СУСПЕНЗИИ В ПЛАСТИНЧАТЫХ (ТРУБЧАТЫХ) ФИЛЬТРАХ
Ключевые слова: Двухфазная суспензия, неизотермическое течение, фильтрование, проницаемость, неньютоновская среда,
механика многофазных сред.
Рассматривается неизотермическое течение двухфазной среды со степенным реологическим уравнением состояния в криволинейных каналах и трубах и процессы фильтрования с образованием осадка в этих областях. Для решения уравнений движения среды применяется метод поверхностей равных расходов. Для определения распределения температуры среды в области течения использовался метод Слезкина. Для вычислении толщины слоя осадка дисперсных включений и их средней концентрации по поперечному сечению области течения используются соответствующие дифференциальные уравнения, полученные из условия баланса твердой фазы.
Keywords: Two-phase suspension, nonisothermal flow, filtration, permeability, non-Newtonian medium, mechanics of multiphase
media.
Nonisothermal flow of a two-phase medium with a power-law rheological equation of state in curvilinear channels and pipes is considered together with processes of filtration involving cake formation within these regions. For solving the equations of the medium motion, the method of surfaces of equal flowrates is applied. The method of Slezkin is applied for determining temperature distribution of the medium in the flow domain. For calculating the thickness of the cake layer as well as dispersed inclusions and their average concentration across the cross section of the flow domain, appropriate differential equations derived from the balance conditions for solid phase are used.
Введение
В процессах химической технологии, экологии и других смежных отраслях промышленности широко реализуется течение двухфазных сред в каналах и цилиндрических трубах с проницаемыми стенками. Такие течения используются, например, в патронных и динамических фильтрах, фильтрах-сгустителях, трубчатых мембранных элементах. Очень важно при выполнении гидродинамических расчетов таких аппаратов учитывать влияние изменения давления, температуры, концентрации и толщины слоя осадка вдоль канала на характер процесса. Изменение температуры суспензии изменяет вязкостные характеристики среды, что приводит к изменению эффективности процесса фильтрования.
В настоящее время эти факторы учитываются достаточно приближенно. Неизотермическое течение гетерогенных сред с неньютоновской реологией в таких каналах недостаточно изучено. Например, работа [1] посвящена математическому моделированию изотермического течения многофазных сред в областях с проницаемыми стенками, в работах [2, 3] рассмотрены вопросы математического моделирования тонкослойных неизотермических течений двухфазных сред по проницаемым поверхностям на основе методов механики многофазных сред [4].
Целью данной работы является построение математической модели неизотермических течений двухфазных суспензий в трубчатых фильтровальных элементах с учетом процесса фильтрования через стенку рабочей поверхности, и численной схемы их расчета.
Теоретическая часть
Рассматривается ламинарное, осесимметри-ческое, неизотермическое и установившиеся тече-
ние гетерогенной среды с твёрдой фазой в проницаемых каналах и трубах различной геометрической формы. Реологическое уравнение состояния неоднородной среды описывается степенной моделью Оствальда де Виля. Уравнения сохранения записываются в ортогональной системе координат, в которой одна из координатных поверхностей совпадает со стенкой канала, а координатные поверхности составляют семейство нормалей к ней. Схема течения представлена на рис. 1.
Рис. 1 - Схема течения
На вход в трубчатый фильтр подается двухфазная среда с входными температурами фаз 71, ВХ и Т2, ВХ, температура стенки трубы (пористой
поверхности) - Тст. На входе в трубчатый элемент задано давление Рвх , на внешней поверхности стенки поддерживается давление РВЬХ . В дальнейшем для простоты решения будем считать, что Т1, ВХ « Т2, ВХ = ТВХ и температуры фаз мало отличаются друг от друга, т.е. Т « Т2 = Т - температура смеси. При значительном превышении длины над поперечным размером канала, в этой выбранной системе координат, порядок продольной скорости
основного потока будет значительно превышать порядок поперечной скорости. После проведения анализа значимости слагаемых, двумерные упрощенные уравнения сохранения механики многофазных сред [4] в размерных переменных в цилиндрической системе координат, связанной с областью двухфазного неизотермического течения, согласно [1,2], запишутся в виде:
д(г-рх | д(г- Р1 -УО = 0 дг дг
(1)
Р1 -[и
-^12,2 +Р1 ^2 дР
ди
дг
+ V
ди
дР
~ =-«1— + П дг ! дг
Ч
дг
- р12,г + Р1рг =
д(г • р2 - и2) + д(г • Р2 -У;) = 0
дг
дг
(2)
(3)
(4)
ди 2 ди 2 1 дР с^ч
Р2 • | и2 + У2 ^Т I = -«2 ^ + ^12,2 + Р2^2 , (5)
дР
-«2 ^^ + ^12,г +Р2 Рг = 0, (6)
(тт дТ т ВТ 1 1 д ( „ дТ 1
СР Р-| + 1= -дт|г-Я-1Г I, (7)
^ дг дг) г дг ^ дг)
«1( г) + «2(г) = 1, (8)
Уравнения фильтрации, сохранения энергии для течений через тонкий слой образовавшегося осадка и через стенку тонких фильтров (при одномерной фильтрации) соответственно записываются в виде:
дг¥,
дг
т = 0 V" =- К
и, V т
дР
Vт дг
V дТт 1 д[
Ср1 -Р1т ^т • = г -А
дг г дг ^
дгV0 "К дР0 _0 = 0 V" =- —__—
дТ„
1т •"
дг
дг
Ср1 -Рю ^0 1 д
0
Vo дг
дТ0 1 д [ дТ0
-=---1 г • А--
дг г дг I дг
(9) (10) (11) (12)
где т1 =
г дг
ди1 ди1
дг дг
"-1
т(Т, «2) -
динамическая эффективная вязкость жидкой фазы, ^ 12 - вектор силы межфазного взаимодействия; и, V -компоненты вектора скорости; р , Т , А, Ср
плотность, температура, коэффициент теплопроводности и теплоемкость смеси; , ¥г - компоненты
вектора массовой силы; 5 - толщина осадка; Рт, Vm - давление и скорость жидкой фазы при фильтрации через стенку; Кт - коэффициент проницаемости стенки, V т - эквивалентная (кажущаяся)
вязкость жидкости в стенке; Р0, Vo - давление и скорость жидкой фазы при фильтрации через слой осадка; К 0 - коэффициент проницаемости слоя
осадка, V0 - эквивалентная (кажущаяся) вязкость жидкости в слое осадка; «1, «2 - объемные концентрации несущей и твердой фаз соответственно; Р1, Р2 - приведенные плотности; рш , Р1т - "размазанные плотности" несущей фазы при её течении в порах осадка и стенки.
Сила межфазного взаимодействия для гетерогенных систем с вязкой сплошной средой, особенно для систем со сложной реологией, в основном определяется вязким трением на межфазной поверхности. Сила межфазного взаимодействия, действующая на дисперсную частицу со стороны потока неньютоновской жидкости, подчиняющейся степенному реологическому закону, будет зависеть от параметра «, критерия Рейнольдса Яе^ , концентрации «2 и от размеров частицы. Коэффициент сопротивления для степенных жидкостей может быть вычислен с помощью эмпирических формул [4-8].
Коэффициент эффективной вязкости т сильно зависит от концентрации твердых включений и от температуры сплошной среды. В работах [9-11] можно найти различные формулы для эффективной вязкости, справедливых при определенных диапазонах изменения а2 и «. Так как коэффициент " мало зависит от температуры, то для реологического закона Оствальда де Виля можно использовать степенную аппроксимацию зависимости от температуры:
т0(Т) = ! , 5 а пк (13)
= 1 +5 ак-пк т(Т) к
«1
к=1
гдеП = Т -Твх . ТСТ - ТВХ
Граничные условия для задачи (1)-(12) следующие:
при г = 0: ди1/ дг = 0, V1 = 0, ди2/ дг = 0,
дТ / дг = 0 (14.1)
при г = Я2 - 5: Р = Р0, и1 = 0 , «^х = V0, и 2 = 0, V2 = 0, (14.2)
Т0 = Т , А-дТ/дг = А10 • дТ0 /дг . при г = ^2 : Р0 = Рт , Vo = Vm , (14.3)
Т0 = Тт , А1т • дТт / дг = А10 - дТ0 / дг ,
при г = щ : Рт = Рвых , Тт = Тст (14.4)
Граничные условия на входе при г = гн :
и = иы (г), и 2 = и 2н (г), (г) = 0, V2„ (г) = 0, Т (г) = Тех , Р(г) = Рвх .
Здесь и 1Н (г), (г), и 2„ (г), V2„ (г) - заданные распределения скоростей жидкой и твердой фаз на входе в трубчатый элемент; Твх, Рвх - заданные значения температуры и давления среды на входе, Рвых - заданное значение давления за стенкой фильтра. Решение системы уравнений (1)-(12) для
1
граничных условий (14) представляет собой сложную проблему, связанную с большими математическими, вычислительными трудностями. В связи с этим, на основе обоснованных предположений, предлагается упрощенная модель процесса.
Для тонкослойных стенок и осадка вместо уравнений фильтрации через пористую стенку и осадок (9)-(12) можно принять следующие приближенные зависимости:
(15)
V" = ~ ' к-5 - РВЫХ), Тт = Т0 = ТСТ
у 1 2
здесь Vу — скорость фильтрования; к, у — коэффициенты проницаемости и кажущейся "вязкости" эффективной среды, которой является стенка фильтра и осадок.
Граничные условия на границе "область течения -поверхность осадка" примут следующий вид:
Уравнения для поверхностей равного расхода определим из (17). Для этого представим интеграл по одной из формул численного интегрирования и затем продифференцируем полученную разностную формулу по г . Если использовать формулу трапеции, то уравнения для определения линии тока примут вид
оГг
= 0 :
йг1к Ф1 к+1 Г1к+1 — г1к йА1к+1
- +----
й! яД^+1 Ахк+1 йг
йг1к+1 йг йг
к = 1, N — 1 (18)
где Д1к +1 = «1(г1к +1и1к +1 + г1ки1к ). уравнение сохранения импульса (2) для жидкой фазы, на линиях тока к (г) (к = 1, N) , после выполнения соответствующих преобразований запишется в виде системы дифференциальных уравнений:
при г = Я2 —5 : и1 = 0 , a1V1 = Vf, и2 = 0,
Т = Т
СТ
(16)
Уравнения системы (1)-(6) являются нелинейными дифференциальными уравнениями в частных производных, которые в квадратурах не интегрируется. В связи с этим нелинейная система (1)-(6) решается с помощью модифицированного метода поверхностей равных расходов [2], предложенного в [11]. В поле течения для каждой фазы (г = 1 — жидкая фаза, г = 2 — твердая фаза) выбираются поверхности равных расходов (трубки тока) гг к (г) (к = 1, N) , форма которых определяется геометрией сечения канала. Поверхности равных расходов гг 1 (г) и гг N (г) совпадают с поверхностями, ограничивающими область течения. Причем, если область течения среды есть круглая цилиндрическая труба, то линия тока гг N (г) совпадает с поверхностью осадка, а линия гг 1 (г) - с осью симметрии - линией - центром трубы. Задача о развитии течения среды сводится к численному определению положений поверхностей равных расходов, скоростей на них и давления. При отсутствии мас-сообменных процессов изменения расхода сплошной фазы между введенными поверхностями определяются следующим образом
й
г
1к
Ф1к (г) = — $ 2 'К'Щ ' г' и1' йг, к = 2, N, (17) йг
г1к—1
и
йи
1 к
а
1 к
йг
1 йР 1 —--+ —
Р йг Г к
д_
дг
Р1
ди1
1к
1к
дг
ди1 к п—1
дг
Р1
РП7, к + , к , к = 1, N — 1 (19)
Аналогичные уравнения получим для второй фазы:
йг21 0 йг2 к+1 йг2 к Ф 2 к+1 г2 к+1 — г2 к йА 2 к+1
—¡— = 0,-=-+----,
йг йг йг жА 2 к+1 А 2 к+1 йг
к = 1, N — 1 (20)
где А2 к+1 = «2 (г2 к+1и2 к+1 + г2 ки2 к ) .
йи 2 к
и
2 к '
йг
а2 йР 1
--.--1--
Р2 йг Р2
Риг, к + , к (21)
к = 1, N — 1 .
Значения поперечных скоростей определяются из соотношений, определяющих линии тока:
йггк
Vк = игк^к
йг
г = 1 и 2, к = 2, N — 1.
При решении системы уравнений (18) - (21) необходимо учитывать влияние температуры на вязкость жидкой фазы[12] . Для этого решение уравнения (7) с целью получения аналитических зависимостей ищется методом Слезкина [13].
Согласно этому методу конвективный член (левую часть) уравнения (7) заменим его средним значением по высоте области течения при помощи функции:
После их дифференцирования по г не трудно показать, что при отсутствии массообмена имеют место соотношения:
Ф1к (г) = 0, к = 2, N — 1,
Ф1N (г) = 2 • ж • «1 • (Я2 — 5) • V (г, Я2 — 5).
В( г) =
с „ 'Р 1 К2—5 ( дТ дТ Л Р И 1 $ I и ■ — + V •д7\йг (22)
2 Я2 —5
0
дг
дг
Подставив в левую часть уравнения (7), получим уравнение:
Проинтегрировав его с уче-
В( ) 1 д ( дТ
В(г) =---1 г--
г дг ^ дг
том граничных условий (14.1) и (16), получим для температуры выражение:
г 2 -(( -5)2 Т (г, г) = Тст + В( г)--^-(23)
В полученном выражении и при решении уравнений течения и фильтрации среды имеются две неизвестные величины В(г) и 5(г, /). Для их определения используем соотношение (22) и уравнение баланса осадка, которое в дальнейшем будет выведено. Из (23) найдем частные производные.
- = В(г)•г , дг 2
дТ = г2 -( -5)2
■В' +
5'-( -5)
•В + Т^Т
дг 4 2
Умножив их на и и V соответственно, сложив и проинтегрировав по г от 0 до ^ -5 получим:
дТ дТ г и • — + V- — = и---В(г) + дг дг 2
(
(24.1)
+ V-
г 2 -((?2 -5)2 -В '+5'( -5
4
2
-В + Т'с т
= V---(42 5) • В' + [у-5 •( 5) + и• В + V• Т'СТ .
Проинтегрируем (24.1) от 0 до Я2 -5, предварительно разбив область интегрирования на подобласти [, гк+1 ] по г . Далее подставив полученное выражение в (22) и выполнив соответствующие преобразования, получим линейное обыкновенное дифференциальное уравнение для В(г):
а-В' + Р- В = у, (24.2)
здесь:
«= гИ^^ + ¿Л^^+11а* ,
^ Ук-Тк + и к+1 •гк+1 + 5'-((?2 -5)-(( + Vk+0
2
2
• Агк -
А (я - 5) N-1
-1 Г = -5(к + Vk+l)Агк-ТСт . (24.3)
с р • Р
Общим решением уравнения (24.2) является:
В (г) =
[У. е•
«
йг + С
Р1«йг
(25)
Для определения константы С в формуле (25) используем условие Т (0, г) = Твх, где Твх есть среднее значение распределения температуры во входном сечении:
В(0) =
^У е ¡р/«йг
йг
+ С
г=0
С = В(0)-е
|Р!«-йг\
е\р!«йг
«
йг (26)
г=0
Т (0, г) = Тст + В(0)
г 2 -((2 - 5(0))2
Так как
5(0) = 0, то: Т (0, г) = Тст + В(0) •
г 2 - Я 2
Найдем среднее значение распределения температуры Т (0, г) в входном сечении:
Твх=—•• |Т(0,г) йг=—- \
Я
2 0
Я
Тст+В(0)
2 021 г - Я
1 1
Я2
В(0)-4
Тст -Я2 +
В(0)
20
( 3 г
3
(и 3
йг =
Я
\\
3
2 Я23
-Тст -В(0)-
Здесь Твх - температура среды на входе в фильтр.
Отсюда: В(0) =
6 • (тст - Твх)
я2
(27)
Подставив (27) в (26), получим:
с _ б- (тст - Твх ) \р1«йг\г 0 - у \Р/« - йг
«
Я3
йг
г=0
и с учетом этого выражения из (25) получим:
е ^«^.¿г + 6-(ст - Твх ) х
■I
В( г) =
Я-
х е
| Р« • йА
у \Р/«• йг
'-/ — •е «
•йг
г=0
е!р/«йг йг + 6- (тст- твх ) • е1р«-йг\ г=
Я3
-! Р1«йг =
-! Р«-йг =
= е-!Р«й •!у.е!Р'«й-йг +6-(тст -Твх\ 0 « Я 3
-! Р« • йг
« • <
х е 0 . Тогда окончательно получится зависимость:
г
(г) = е ' - ! —-е1
В(г) = е~!Р«йг -!У- е!Р«йг йг + 0 «
6-(ст - Твх) е"! Р«йг
Я23
(28)
Во всех полученных выражениях пока имеется неизвестная величина 5(г, t).
Выведем уравнение для вычисления накопления осадка 5 = 5(г, t), меняющейся как по длине канала, так и по времени. Для этого построим уравнение баланса для толщины слоя осадка твердой фазы в элементарном объеме осадка, ограниченном сечениями г иг + Аг , и двумя поверхностями, полученными от вращения линий г = Я2 и г = Я2 -5(г) вокруг оси фильтра. В моменты времени t и t + Аt соответствующие элементарные объемы вычисляются по формулам:
4
4
4
Я
0
Я
4
6
е
Р
е
□(г + Аг) = 2ж' Я2 А 5(г, г + Аг) и П(г) = 2ж' Я2 А 5(г, г).
Изменение количества твердых частиц в рассматриваемом объеме определяется потоком оседающих частиц на поверхность осадка за
время Аг (движением осажденной массы по направлению оси 2 пренебрегаем).
а20 (( + Аг) — □(/)) = . (29)
Площадь поверхности осадка в элементарном объеме это площадь поверхности вращения кривой Я2 —5(г, г) вокруг оси 02. Эта площадь определится формулой: Б5 = 2ж' (Я2 — 5(г, г))' лД + 5'22 ' Аг . Объем осадка,
накопившегося на этой поверхности за время Аг определится формулой:
Оп = 2ж'(Я2 —5)'л/1 + 52Аг«2 х xV2( г, Я2 —5)' Аг (30.1)
П(г + Аг) — П(г) = 2ж - Я2' Аг' (5(г,г + Аг) — 5(г,г)). (30.2) Тогда из (30.1) и (30.2) получим соотношение: 2ж'Я2 «20 'Аг'(5(г,г + Аг) — 5(г,г)) =
= 2ж'(Я2 — 5)^1 + 5'/ Аг«2 ^2(г,Я2 —5)' Аг
Сократив последнее равенство на Аг , далее разделив на «20 ' 2ж' Я2 ' Аг, получим:
5( г,г + Аг) — 5( г,г) = «2'(Я2 —5) Г5 ^ я2 —5)
Аг
«20 'Я2
Выполнив предельный переход при Аг ^ 0 , получим уравнение баланса массы для слоя осадка
^ = «2 '(Я2 —5)' /1+5,2^я2 —5) (31) дг а20 'Я2 Для средней по поперечному сечению канала объемной концентрации твердой фазы, присутствующей в построенной математической модели, из балансовых соотношений получим уравнение: йа2
йг
= 2ж' (( — (г, Я2 —5) «2/ [22 (гм) +
+ $ 2ж' (Я2 (г, Я2 —5)' йг
(32)
Интегральное условие сохранения сплошной фазы для фильтра длиной Ь запишется в виде:
Я2 5
$ 2ж'« ' и1(г, г)' г ' йг +
+ $2ж'(Я2 —5)^/ (г, Я2 —5)' йг =
(33)
Уравнение (33) позволяет определить длину фильтровального элемента, необходимой для полного разделения исходной суспензии при условии равенства нулю первого слагаемого. Второе слагаемое уравнения (33) определяет количество фильтрата при заданной длине фильтровального элемента Ь.
Математическая модель задачи представляет собой систему, состоящей из уравнений (13), (15), (18)-(21), (23), (28), (31) -(33) , которые должны решаться совместно.
Вычислительный алгоритм
1. В начальный момент времени г = 0 осадок отсутствует 5(г) = 0 .
2. Во входном сечении г0 = гн : задаются значения давления Рвх , температуры суспензии Твх , значение эффективной вязкости ш0) (Твх ), выбираются распределения положений поверхностей равных расходов гг к (г) (к = 1, N) и соответствующие значения скоростей игк (г = 1,2).
3. Для следующего сечения г у+1 = г у + Аг
численным методом равных расходов решается система обыкновенных дифференциальных уравнений (18)-(21), совместно с (13), (15), (23), (28). При этом для сечения г = г у, у = 0,1,2,... по (18)-(21) вычисляются значения правых частей этих дифференциальных уравнений с использованием апроксима-ций выражений для скоростей по полной системе координатных функций с учетом граничных условий. Значение правой части дифференциального уравнения для давления, входящего в эту систему, определяется методом прогонки . В результате находятся искомые функции гг к (г) (к = 1, N) и иг к (г = 1,2), давление Р(г), температура Т (г, г^), эффективная вязкость т(Т)
на следующем шаге интегрирования по продольной координате. По найденному значению давления, из системы уравнений (15), (32) вычисляются скорость оттока сплошной фазы Vу, а 2 .
4. Процедура п. 3 повторяется до достижения конечного сечения канала г = Ь .
5. Из уравнений (31) определяются распределения толщины осадков вдоль канала для следующего момента времени г = г + Аг < гк.
6. Для нового момента времени повторяются пункты 2-5.
Заключение
Проведена модификация метода равных расходов для расчета процессов течения и фильтрования двухфазной среды в трубчатых фильтрах в неизотермических условиях. Построенная математическая модель позволяет рассчитать процессы неизотермического течения двухфазных сред в каналах и трубах с проницаемыми стенками и определить все характеристики процесса разделения и размеры фильтра.
Литература
1. Ибятов Р.И., Холпанов Л.П., Ахмадиев Ф.Г., Бекбула-тов И.Г. ТОХТ., 41, 5, 514-523, (2007).
2. Ахмадиев Ф.Г., Фазылзянов Р.Р., Галимов Р.А. ТОХТ., 46, 6, 620-630, (2012).
3. Ахмадиев Ф.Г., Фазылзянов Р.Р., Галимов Р.А. Вестник Казанского Технологического Университета, 16, 14, 101-105, (2013).
4. Нигматулин Р.И. Динамика многофазных сред. Ч.1. М., 1987.
г
н
г
н
5. Келбалиев Г.И. ТОХТ., 45, 3, 264, (2011).
6. Потанин А.А., Черномаз В.Е., Тараканов В.М., Урьев Н.Б. ИФЖ., 60, 1, 32, (1991).
7. Acharya A., Mashelkar R.A., Ulbrecht J. Reol. Acta., 15, 9, 454, (1976).
8. Tanaka H., White J.L. J. of Non - Newt. Fluid. Mech., 7, 4, 333,( 1980).
9. Шмаков Ю.И., Шмакова Л.М. В кн.: Механика жидкостей и газа. Ташкент, 1980. С.77.
10. Иванов В.А. В кн.: Процессы тепло- и массопереноса вязкой жидкости. Уральский научный центр АН СССР. 1986. С.79.
11. Холпанов Л.П., Шкадов В.Я. Гидродинамика и тепломассообмен с поверхностью раздела. М.: Наука, 1990.
12. Кадыйров А.И., Вачагина Е.К. Вестник Казанского Технологического Университета, 16, 19, 49-53, (2012).
13. Тарг С.М., Основные задачи теории ламинарных течений. - М.: ГИТТЛ, 1951. - 420 с.
© Ч. Х. Исянов - асп. каф. процессов и аппаратов химической технологии КНИТУ, инженер ООО ИВЦ «Инжехим», [email protected]; Ф. Г. Ахмадиев - д-р техн. наук, проф., зав. каф. прикладной математики, КГАСУ; М. И. Фарахов- д-р техн. наук, проф. каф. процессов и аппаратов химической технологии КНИТУ, дир. ООО ИВЦ «Инжехим», [email protected]; И. Г. Бекбулатов - доцент кафедры прикладной математики КГАСУ.
© C. H. Isyanov - graduate student of the Department of Processes and Apparatuses of Chemical Engineering of KNRTU, engineer of LLC EPC «Inzhekim», [email protected]; F. G. Ahmadiev - Head of the Department of Applied mathematics, Doctor of Technical Sciences, Professor of Kazan State University of Architecture and Engineering; M. I. Farakhov - Doctor of Technical Sciences, Professor of the Department of Processes and Apparatuses of Chemical Engineering of KNRTU, Director of LLC EPC «Inzhekim», [email protected]; I. G. Bikbulatov - Associate Professor, Department of Applied Mathematics of Kazan State University of Architecture and Engineering.