Mathematical modelling of mass-transfer under the evaporation of multicomponent liquids

Full Text

Процесс испарения многокомпонентных жидкостей зависит от многих факторов: состава жидкой фазы, изменения ее объёма, наличия градиентов температур и концентраций, взаимного влияния процессов теплообмена и массообмена и т.д. Многообразие факторов, влияющих на процесс, обуславливает многообразие подходов к данной проблеме [1–6]. Однако среди значительного количества работ практически отсутствуют исследования процесса испарения многокомпонентных смесей для открытых систем с учётом массопереноса веществ в жидкой фазе, особенно если жидкость содержит неиспаряющиеся (или плохо испаряющиеся) компоненты с высокой концентрацией. Впервые эта проблема исследована в работе [1], где рассматривается испарение бинарной смеси, один из компонентов которой не испаряется. Решение для смеси, содержащей N компонентов, с учётом теплообменных процессов приведено в [2], но без учёта наличия в растворе неиспаряющихся веществ. В работе [3] дано полное математическое описание испарения в присутствии неиспаряющегося компонента с учётом движущейся во времени границы. Следует отметить, что в указанных работах [1–3] испарение рассматривается с учётом диффузионного массопереноса в газовой и 98 Математическое моделирование массопереноса в процессе испарения . . . жидкой фазах при отсутствии конвективных процессов в газовой фазе. При движении воздушных масс над поверхностью зеркала испарения применимость данных работ будет весьма ограничена, в то же время существующие на сегодня модели, отражающие данную специфику, не учитывают массопереноса в жидкой фазе [4–6]. В связи с этим в задачу настоящей работы входило изучение закономерностей процесса испарения многокомпонентных жидкостей в открытых системах при постоянной температуре и скорости движения газовой фазы. Особенность данного исследования заключается в разработке математической модели процесса одномерного испарения многокомпонентных жидких смесей с высокой концентрацией одного из компонентов (далее по тексту — растворитель) и с низкой концентрацией всех остальных N компонентов (далее по тексту — растворённые вещества) при постоянном движении воздушного потока над границей раздела фаз. Будем рассматривать два случая: когда растворитель испаряется и не испаряется. При разработке математической модели будем исходить из следующих допущений: – в начальный момент времени концентрации веществ распределены равномерно; – испарение веществ с поверхности жидкости подчиняется закону Герца— Кнудсена; – парциальное давление насыщенного пара каждого компонента над раствором подчиняется закону Рауля; – массоперенос растворённых веществ в жидкой фазе осуществляется за счет диффузионных процессов, подчиняющихся законам Фика при постоянном коэффициенте диффузии; – массоперенос растворителя в жидкости отсутствует; – интенсивность испарения растворителя в начальный момент времени меньше интенсивностей испарения каждого компонента из N растворенных веществ. При этом будем учитывать уменьшение объёма жидкости за счёт испарения путём определения скорости движения границы, разделяющей жидкость и газ. Математическое описание испарения многокомпонентной жидкости при постоянной температуре растворённых веществ и отсутствии испарения растворителя имеет следующий вид: ∂2 ∂ Ci (x, t) = Di 2 Ci (x, t) при h(t) < x < ∞; ∂t ∂x Ci (x, t) = C i = const при t = 0; ∂ d Ci (x, t) h(t) + Di Ci (x, t) = Ji (t) при x = h(t); dt ∂x Ci (x, t) = C i = const при x = ∞; d h(t) = dt N j=1 1 Jj (t) ; ρj h(t) = 0 при t = 0; (1) (2) (3) (4) (5) (6) 99 А. Ю. Ш у п л я к, С. П. Ш к а р у п п а, А. М. Ш т е р е н б е р г Jj (t) = αYj (t)pj N Yj (t) = µj µj /T ; (7) Ci (x, t) при x = h(t), (8) −1 Ck (x, t)/µk k=0 где t — время, с; x — координата в жидкой фазе, перпендикулярная границе раздела фаз (начало координат соответствует нахождению границы раздела фаз в начальный момент времени, см. рис. 1), м; i — порядковый номер компонента (i = 0 соответствует растворителю), i = 1, 2, . . . , N ; C i , Ci — начальная концентрация и концентрация в произвольный момент времени (г/м3 ) i-того компонента; Di , Ji , µi , pi , ρi — коэффициент диффузии (м2 /c), интенсивность испарения (г/(см2 )), молярная масса (г/моль), давление насыщенных паров (Па) и плотность (г/м3 ) i-того компонента соответРис. 1. Схема движения границы: 1 – слой испаряющейся жидко- ственно; Yi — мольная доля компонента в рассти; 2 – слой жидкости, который ис- творе на границе раздела фаз; h(t) — функпарится за время t ция зависимости координаты границы раздела фаз от времени, м; N — число растворённых веществ в растворителе; α — эмпирический коэффициент, зависящий от движения воздушных масс над поверхностью жидкости. В задаче (1)–(8) соотношение (2) представляет начальное условие, которым задаётся равномерное распределение концентрации всех компонентов в растворе в момент времени t = 0. Соотношение (3) представляет граничное условие, согласно которому на границе раздела фаз сумма диффузионного потока, определяемого первым законом Фика [7], и потока, возникающего за счёт движения границы, i-того испаряющегося растворённого вещества равна его интенсивности испарения. Соотношение (4) характеризует равенство концентрации i-того компонента начальному значению при x → ∞. Величина dh(t)/dt — скорость конвективного потока (стефановский поток [7]), возникающего от поступления веществ к поверхности раздела фаз за счёт изменения объёма жидкости. Выражение (6) задаёт координату границы раздела фаз «газ–жидкость» в начальный момент времени. Выражение (7) представляет собой совмещение законов Герца—Кнудсена и закона Рауля, где Yi pi — парциальное давление насыщенного пара испаряющегося компонента над раствором. Краевая задача (1)–(8) представляет систему N дифференциальных уравнений с нелинейным граничным условием третьего рода (совокупность выражений (3), (7) и (8)) и соотношениями (5)–(6), описывающими движение верхней границы. Для упрощения процесса получения аналитического решения задачи (1)– (8) перейдём к подвижной системе координат и приведём граничное условие (3), (7), (8) к линейному виду. Для перехода к подвижной системе координат с фиксированной верхней границей введём новую независимую переменную z = x − h(t) [9]. Считая, что концентрации растворённых веществ малы по 100 Математическое моделирование массопереноса в процессе испарения . . . сравнению с концентрацией растворителя, соотношение (7) в новой системе координат приведём к виду Ji (t) = ki Ci (z, t) при z = 0, (9) где ki — константа скорости (м/с) испарения компонента, определяемая по формуле ki = αµ0 pi / C 0 µi T ; (10) µ0 , C 0 — молярная масса (г/моль) и начальная концентрация (г/м3 ) растворителя. Введём величину b как среднее значение стефановского потока за некоторый промежуток времени tmax : b= tmax 1 tmax 0 d 1 h(t)dt = dt tmax N i=1 1 ρi tmax Ji (t)dt . (11) 0 С учётом введения подвижной системы координат и выражений (9) и (11), систему уравнений (1)–(8) можно представить в виде ∂ ∂2 ∂ Ci (z, t) = Di 2 Ci (z, t) + b Ci (z, t) при 0 < z < ∞, ∂t ∂z ∂z Ci (z, t) = C i = const при t = 0, ∂ Di Ci (z, t) = (ki − b)Ci (z, t) при z = 0, ∂z Ci (z, t) = C i = const при z = ∞, (12) (13) (14) (15) где i = 1, 2, . . . , N , т.е. в виде N независимых дифференциальных уравнений (12) с граничными условиями (13)–(15) для каждого растворённого вещества. Решение дифференциального уравнения (12) при заданных граничных условиях в интегральной форме приведено в [9]: ∞ Ci (z, t) = C i Gi (z, ψ, t)dψ, 0 где Gi (z, ψ, t) — функция Грина, определяемая выражением b(ψ − z) b2 t (z − ψ)2 1 Gi (z, ψ, t) = √ exp − exp − + 2Di 4Di 4Di t 2 πDi t (z + ψ)2 2ki − b ∞ (z + ψ + η)2 2ki − b + exp − − exp − − η dη . 4Di t Di 4Di t 2Di 0 После преобразований решение дифференциального уравнения (12) примет вид Ci (z, t) = Ci z + bt ki − b bz z − bt erfc − √ − exp − erfc √ + 2 ki Di 2 Di t 2 Di t 101 А. Ю. Ш у п л я к, С. П. Ш к а р у п п а, А. М. Ш т е р е н б е р г + z + (2ki − b)t 2ki − b (ki − b)(ki t + z) √ exp erfc ki Di 2 Di t , 0 t tmax . (16) Подстановкой выражения (16) при z = 0 в выражение (9) находим функцию интенсивности испарения от времени при условии, что стефановский поток постоянен (b = const): Ji (t) = C i b t b + erfc − 2 2 Di 2ki − b (ki −b)ki t/Di + erfc e 2 ki − b 2 t Di , 0 t tmax . (17) Следуя [7], преобразуем второе слагаемое выражения (17) для вычисления при больших значениях времени: применим разложение функции erfc(x) в асимптотический ряд. После упрощений получим выражение Ji (t) = C i b b erfc − 2 2 t + Di b2 Di − 4Dt i e πt ∞ n 2n Di (2ki − b)2n tn n=0 n (1 − 2j) . j=0 При больших значениях времени t имеем [7] ∞ n=0 n 2n Di (2ki − b)2n tn erfc(−0,5b t/Di ) ≈ 2, n (1 − 2j) ≈ 1, (18) j=0 exp(−0,25b2 t/Di ) ≈ 1. (19) С учётом выражений (18), (19) получаем приближённую формулу для расчёта интенсивности испарения при больших значениях времени t: Ji (t) ≈ bC i + C i Di /(πt). (20) Схематический график зависимости логарифма функции испарения (17) от времени представлен на рис. 2. Анализ результатов, представленных на рис. 2, позволяет заключить, что в начальный момент времени (t = 0) интенсивность испарения растворённых веществ определяется только константой испарения (ki ) (по аналогии с химическими процессами эту область можно назвать кинетической); при больших значениях времени (t → ∞) процесс испарения растворённых веществ можно рассматривать как стационарный: Ji (t) = ki C i при t = 0, (21) Ji∗ , (22) lim Ji (t) = bC i = t→∞ где ki — константа испарения, определяемая по формуле (10); Ji∗ — установившаяся с течением времени интенсивность испарения i-того компонента. Равновесные концентрации растворённых веществ на границе раздела фаз (г/м3 ) при установившемся стационарном режиме определяются подстановкой в формулу (9) выражения (22): Ci∗ = Ji∗ /ki = bC i /ki . 102 (23) Математическое моделирование массопереноса в процессе испарения . . . Выражения (17)–(22) выведены при условии постоянства стефановского потока (формула (11)), но при больших значениях времени осреднения средняя скорость стефановского потока может стремиться к нулю. Рассмотрим принципиальную возможность установления стационарного состояния, для которого характерно при больших значениях времени равенство значений функции её среднему значению за этот же промежуток времени. При подстановке в формулу (11) зависимости (20) получаем следующее соотношение: b= lim tmax →∞ 1 tmax N i=1 1 ρi tmax bC i + C i 0 Di dt πt N = i=1 bC i . ρi (24) Из выражения (24) получаем условие перехода процесса в стационарное состояние: N y i = 1. i=1 Здесь y i = C i /ρi — объёмная доля i-того растворённого вещества на границе раздела фаз в начальный момент времени (t = 0). Сумма объёмных долей испаряющихся веществ при наличии неиспаряющегося испарителя не может быть равна единице, следовательно, предположение о возможности установления стационарного состояния без участия в испарении растворителя было ложным. В этом случае при больших значениях времени осреднения величина b → 0, следовательно, для бесконечного слоя при отсутствии испарения растворителя интенсивность испарения каждого растворённого вещества будет асимптотически стремиться к выражению, характеризующему диффузионную область протекания процесса: Ji (t) ≈ C i Di /(πt). (25) При участии в испарении плохо испаряющегося растворителя его концентрация на границе раздела фаз будет увеличиваться, следовательно, задачу Рис. 2. Схематический график зависимости логарифма интенсивности испарения от времени для слоя бесконечной глубины 103 А. Ю. Ш у п л я к, С. П. Ш к а р у п п а, А. М. Ш т е р е н б е р г (1)–(8) необходимо дополнить выражением для определения объёмной концентрации растворителя через концентрации растворённых веществ: N y0 (x, t) = C0 (x, t)/ρ0 = 1 − yi (x, t), (26) i=1 где yi (x, t) — объёмные доли растворённых веществ (при i = 1, 2, . . . , N ) и растворителя (при i = 0) соответственно в момент времени, отличный от начального при координате x. С учётом испарения растворителя выражения (5) и (11) примут вид d h(t) = dt b= N 1 tmax i=0 N 1 Ji (t) , ρi i=0 (27) tmax 1 ρi Ji (t)dt , (28) 0 где J0 (t) — интенсивность испарения растворителя (г/(с м2 )); ρ0 — плотность растворителя (г/м3 ); отношение J0 (t)/ρ0 — вклад в скорость стефановского потока процесса испарения растворителя (м/с). Рассуждая аналогично предыдущему случаю, с учётом выражений (27) и (28) получаем модифицированное выражение (24) в виде формулы N ∗ C i /ρi + k0 C0 /ρ0 , b=b (29) i=1 ∗ где C0 — равновесная концентрация растворителя на границе раздела фаз 3 ); член k C ∗ /ρ характеризует испарение растворителя в кинетической (г/м 0 0 0 области. ∗ Определяя отношение C0 /ρ0 из выражения (26) и подставляя его в выражение (29), находим N N y i + k0 1 − b=b yi . (30) i=1 i=1 Для нахождения объёмных долей растворённых веществ разделим концентрацию из выражения (23) для каждого компонента на его плотность: yi = by i /ki . (31) Подставляя (31) в (30), после преобразований получаем уравнение относительно величины стефановского потока: N N b 1− y i + k0 y i /ki = k0 . i=1 i=1 Из выражения (32) получаем значение стефановского потока: N b= i=0 104 yi ki −1 . (32) Математическое моделирование массопереноса в процессе испарения . . . При N y0 k0 i=0 yi ki скорость стефановского потока определяется константой испарения растворителя: b ≈ k0 . Для слоя конечной глубины H имеем задачу, аналогичную рассмотренной задаче (1)–(8), где вместо граничного условия (4) необходимо использовать соотношение ∂ C(x, t) = 0 при x = H. ∂x Проведя аналогичные операции, получим задачу N независимых дифференциальных уравнений (12) с граничными условиями (13)–(15) для каждого растворённого вещества, где вместо граничного условия (15) принимается соотношение ∂ C(z, t) = 0 при z = H − h(t). (33) ∂z Из граничного условия (33) следует, что для ограниченной области жидкой фазы при подвижной системе координат необходимо учитывать движение нижней границы. Между тем для начальной стадии процесса испарения слой можно считать бесконечным. С течением времени концентрации растворённых веществ на поверхности раздела фаз снижаются: уменьшается скорость стефановского потока. При достижении изменения концентрации i-того компонента влияние движения границы может оказаться незначительным. Введём параметр эффективной толщины слоя N tmax 1 ρi Heff = H − i=0 Ji (t)dt . (34) 0 С учётом выражения (34) граничное условие (33) приводится к соотношению ∂ C(z, t) = 0 при z = Heff . ∂z Таким образом, с помощью упрощающего предположения (34) приводим N дифференциальных уравнений (12) с подвижной нижней границей к классическому виду. Решение дифференциального уравнения (12) при заданных граничных условиях в интегральной форме приведено в [9]: Heff Ci (z, t) = C i Gi (z, ψ, t)dψ, 0 где Gi (z, ψ, t) — функция Грина, определяемая так: Gi (z, ψ, t) = exp b(ψ − z) b2 t − 2Di 4Di ∞ n=0 1 2 fn (z)fn (ψ) exp(−Di θn t); Bn fn (z) = cos(θz) + (2ki − b) sin(θz)/(2Di θn ); 105 А. Ю. Ш у п л я к, С. П. Ш к а р у п п а, А. М. Ш т е р е н б е р г Bn = − 2 2 b 4Di θn + (2ki − b)2 2ki − b Heff Heff (2ki − b)2 + ; + + 2 2 2 2 2 2 4Di θn 4Di θn 2 4Di θn + b2 8Di θn θn — множество положительных корней трансцендентного уравнения 4Di (ki − b) tg(Heff θn ) = 2 θ 2 + b(2k − b) . θn 4Di n i После преобразований решение дифференциального уравнения (12) примет вид zλ −H b eff ∞ ξn ξn cos Ci (z, t) = 2C i e n=0 zξn Heff + λΣ sin zξn Heff 2 2 2 e−Di t(ξn +λb )/Heff 2 2 2 ξn − λb λΣ + (ξn + λ2 )(ξn + λ2 )/(λΣ − λb ) Σ b , (35) где λΣ , λb , λk — критериальные числа, определяемые выражениями λΣ = λk − λb ; λb = Pe/2; λk = ki Heff /Di ; (36) Pe = bHeff /Di — диффузионный критерий Пекле; ξn — множество положительных корней трансцендентного уравнения 2 tg(ξn ) = ξn (λΣ − λb )/(ξn + λb λΣ ) (37) на отрезках [nπ; nπ + π/2], n = 0, 1, 2, . . . . Подставляя значение концентрации вещества в жидкости на границе раздела фаз, определяемой из выражения (35) при z = 0, в выражение (9), после преобразований получаем выражение для расчёта интенсивности испарения: J(t) = 2 Di C i Heff ∞ 2 2 2 λk ξn exp − Di t[ξn + λ2 ]/Heff b . 2 2 ξ 2 − λb λΣ + (ξn + λ2 )(ξn + λ2 )/(λΣ − λb ) Σ b n=0 n (38) С течением времени ряд (38) сходится к первому члену (n = 0) и функция интенсивности испарения приводится к экспоненте (см. рис. 3): Ji (t) = A1 exp(−A2 t), где A1 = 2 1 2CDi λk ξ0 (λΣ − λb ) , 2 2 2 Heff (λΣ − λb )(ξ0 − λb λΣ ) + (ξ0 + λ2 )(ξ0 + λ2 ) Σ b 2 2 A2 = Di (ξ0 + λ2 )/Heff . b Для того чтобы математическая модель, выраженная формулами (37) и (38), имела смысл, необходимо, чтобы у трансцендентного уравнения (37) на отрезке [0, π/2] был корень. Запишем условие нахождения корня на данном отрезке: (λΣ − λb )/λb λΣ > 1. 106 Математическое моделирование массопереноса в процессе испарения . . . Рис. 3. Зависимость логарифма интенсивности испарения от времени для слоя глубины H Подставляем в это условие вместо λΣ выражение (36) и решаем неравенство относительно λb с учётом, что λk и λb больше нуля и λb < λk : λb < 1 + 0,5λk − 0,5 λ2 + 4. k Последнее выражение можно переписать в виде Pe < 2 + λk − В предельном случае lim [λk − λk →∞ λ2 + 4. k (39) λ2 + 4] = 0 и Pe < 2. k При невыполнении условия (39) процесс не выходит на экспоненту, изменение интенсивности испарения вещества во времени будет описываться выражением (17), как для слоя жидкости бесконечной глубины. Если растворитель не испаряется или значение числа Пекле настолько мало, что им можно пренебречь (λb → 0), то выражение для расчёта интенсивности испарения примет вид Ji (t) = 2λ2 k Di C i Heff ∞ 2 2 exp(Di ζn t/Heff ) , 2 λ2 + λk + ζn k n=0 (40) где ζ определяется из решения трансцендентного уравнения tg(ζ) = λk /ζ. (41) Из формул (17), (22), (25), (37), (38), (39), (40) и (41) для расчёта испарения следует, что если растворитель не испаряется, то конвективный (стефановский) поток, возникающий за счёт поступления веществ к границе раздела фаз, с течением времени прекращается, и испарение будет стремиться 107 А. Ю. Ш у п л я к, С. П. Ш к а р у п п а, А. М. Ш т е р е н б е р г к диффузионному режиму (выражения (40) и (41), при этом в начальные моменты времени, пока концевой эффект слоя не проявляется, для расчёта можно использовать выражение (25)). При испарении растворителя для веществ, диффузионный критерий которых соответствует условию (39), может установиться квазистационарное состояние, выраженное формулой (22). Длительность стационарного состояния обуславливается соблюдением условия (39). Когда условие (39) нарушится, для некоторого растворённого вещества его испарение будет описываться формулами (37) и (38), а процесс будет стремиться к стационарному испарению без участия этого компонента. Выводы. Таким образом, на основании существующих модельных представлений получена математическая модель, описывающая испарение многокомпонентной жидкости, содержащей в избытке один из компонентов (растворитель), интенсивность испарения которого меньше интенсивности испарения остальных компонентов (растворённых веществ). Исходя из модели при рассмотрении процесса испарения многокомпонентных жидкостей можно выделить три группы веществ: растворитель, растворённые вещества, для которых выполняется условие (39), и растворённые вещества, для которых условие (39) не выполняется. Данная модель характерна для смесей, содержащих значительное количество плохо испаряющихся компонентов примерно одинаковой плотности или разной плотности, но высокой вязкости. С течением времени концентрация быстро испаряющихся компонентов на границе раздела фаз падает, а плохо испаряющиеся компоненты накапливаются. Таким образом, даже если в начальный момент времени концентрации легко испаряющихся компонентов высоки, с помощью разработанной математической модели можно рассчитать конечную стадию испарения многокомпонентных открытых систем. Введённая величина b может рассматриваться как конвективный поток, возникающий при воздействии на жидкость сторонних сил, а стефановский поток — как вклад в эту величину. В этом случае при отсутствии испарения растворителя система также может достигнуть квазистационарного состояния, но при этом следует учитывать массоперенос растворителя. Полученная математическая модель процесса испарения многокомпонентных жидкостей позволяет прогнозировать состояние открытых систем в различные интервалы времени при изменении концентраций веществ в широких диапазонах.

About the authors

Anatoliy Yu Shuplyak

Samara State Technical University

Email: yuashuplyak@mail.ru
Engineer, Dept. of Chemical Engineering & Industrial Ecology. 244, Molodogvardeyskaya st., Samara, 443100, Russia

Svetlana P Shkaruppa

Samara State Technical University

Email: ecology@samgtu.ru
(Ph. D. Chem.), Associate Professor, Dept. of Chemical Engineering & Industrial Ecology. 244, Molodogvardeyskaya st., Samara, 443100, Russia

Alexander M Shterenberg

Samara State Technical University

Email: ashter53@mail.ru
(Dr. Phys. & Math. Sci.), Professor, Head of Dept., Dept. of General Physics and Physics of Oil and Gas Production. 244, Molodogvardeyskaya st., Samara, 443100, Russia

References

Supplementary files