INVESTIGATION OF STATIC AND DYNAMIC REGIMES OF MASS TRANSFER PROCESSES

Аннотация

A mathematical model for mass transfer process in rectification columns of poppet type is obtained. At the result of the investigation the ordinary differential equations system is obtained with the transition to equations in partial deriva- tives. Investigations of processes of detection establishment of mass transfer processes in dynamic regimes are per- formed. Static and dynamic characteristics for industrial rectification columns are calculated.

Толық мәтін

Для математического описания процессов массо- обмена в ректификационных колоннах широко ис- пользуются системы уравнений в частных производ- ных [1; 2]. Такое описание вполне естественно для колонн насадочного типа, но требует отдельного обоснования для тарельчатых колонн, так как в по- следнем случае объект по своей природе дискретен. В [2] был развит подход, основанный на детальном рассмотрении процессов для отдельной тарелки. На основе физических представлений о гидродинамике жидкости в тарелке и барботаже парового потока бы- ли получены уравнения баланса массы с учетом фазо- вого перехода компонентов [2]. В результате прове- денных исследований построена система обыкновен- ных дифференциальных уравнений, а затем и в част- ных производных, с использованием формул разло- жения Тейлора в предположении малого различия параметров потока на соседних тарелках. Полученная система уравнений в частных производных для та- рельчатых колонн имеет такой же вид, что и для наса- дочных, благодаря чему возможно описание колонн различного типа с помощью единого математического аппарата и сравнительно несложного пересчета пара- метров для колонн различных типов. 36 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева Рассмотрим систему уравнений, описывающую где un , vn – собственные функции частот; λn – ком- массообмен для бинарной смеси в ректификационной колонне тарельчатого типа [2]: плексные числа; cn – произвольные коэффициенты. Для применения конечно-разностной методики ∂ ((HS + S h ( H ))nx ) − Δl ∂ (Qnx ) = к решению системы (1) введем следующие обозначе- ∂t 1 2 k πD3/ n ∂l ( y − kx ) + Φ1 , ния: u = (HS + S h ( H ))nx , v = , k y = k πD3n . 6τ 6τ ∂ ∂ ⎛ k πD3/ n y ⎞ (1) Линеаризуя (1) в окрестности стационарных парамет- ров, получим систему без учета внешнего воздейст- (S2 (Δl − h0 ) ny ) + Δl т = ⎜ ⎟ вия: ∂t ∂l ⎝ 6τ ⎠ ∂u ∂u k πD3/ n (kx − y ) + Φ2 , * ∂t ∂l (3) 6τ ∂v + B ∂v = D (k *u − v ), где x, y – массовые, а n , n – мольные концентрации компонентов в паре и жидкости соответственно; H, h(H) – уровень жидкости в переливном патрубке и ∂t ∂l где A, B – величины, пропорциональные скоростям потоков. Граничные условия с учетом рециркуляции на тарелке; Q – поток жидкой фазы; Δl – расстояние взаимодействующих для этой системы можно запи- между тарелками; S1 , S2 , kт , D – геометрические сать в виде [2]: параметры тарелки; τ – величина, характеризующая α du + β dv = γ u + δ v, l = 0, время образования пузырька пара в жидкости; 1 dt 1 dt 1 1 (4) Φ1 , Φ2 – функции внешнего воздействия, обуслов- α du +β dv = γ u + δ v, l = L, ленные вводом и выводом потоков сырья и целевых 2 dt 2 dt 2 2 продуктов. Большой интерес представляет исследование раз- где αi , βi , γi , δi (i = 1, 2 ) – коэффициенты, зависящие личных нестационарных режимов, в частности про- цессов установления. При малых отклонениях от рав- новесного состояния возможно проведение линеари- зации системы и использование метода стоячих волн от параметров потоков в кубе (l = 0) и дефлегматоре (l = L). Конкретный вид этих коэффициентов приве- ден в работе [2]. Решение задачи (3), (4) ищем в виде для исследования спектра собственных частот. Не- смотря на простоту, этот метод позволяет не только u = u (l ) eλt , v = v (l ) eλt . получить качественное представление о характере процесса установления, но и определить важные ко- личественные характеристики времени установления, резонансные свойства системы, а также выявить об- Подставим эти решения в систему (3) и граничные условия (4): u (l ) λ− Au′ (l ) = D (v (l ) − k *u (l )), (5) ласти неустойчивости в пространстве параметров сис- темы. v (l ) λ+ Bv′ (l ) = D2 (k *u (l ) − v (l )), Вычисление времени установления необходимо по двум причинам: во-первых, это важно для прогнози- рования времени перехода с одного стационарного режима к другому при изменении скорости поступле- α1λu (0 ) + β1λv (0 ) = γ1u (0 ) + δ1v (0 ), α2 λu (0 ) + β2 λv (0 ) = γ2 u (0 ) + δ2 v (0 ), Решим систему (5): l = 0, l = L. (6) ния сырья или его состава; во-вторых, это нужно для u (l ) = aeμl , v (l ) = beμl . оптимального проведения расчетов в более сложных программах решения системы (1) по конечно- разностной методике. В этом случае программа вы- числения собственных частот и времени установления включается как блок в общую программу и позволяет Подставляя u(l) и v(l) в (3), получим систему на собст- венные значения λ: aλ− Aaμ = D (b − k *a ), выбрать оптимальный шаг интегрирования по време- ни, совместимый с устойчивостью и удовлетвори- bλ+ Bbμ = D2 (k *a − b), тельной аппроксимацией. Для контроля точности рас- из которой найдем выражения для μ1,2 (λ ) : четов по конечно-разностной методике в качестве 1 ⎛ λ λ D D k * ⎞ тестовых также могут быть использованы аналитиче- μ1,2 = ⎜ − − 2 + 1 ⎟ ± ские решения, представленные в виде суперпозиции 2 ⎝ A B B A ⎠ стоячих волн. Если известно N собственных значений, то можно 1 ⎛ λ λ D D k * ± − − 2 + 1 ⎞ λ2 + λD λk * D + 2 + 1 . ⎜ ⎟ выписать соответствующее N параметрическое се- мейство решений в виде 4 ⎝ A B B A ⎠ AB AB AB ⎛ u ⎞ N = c ⎛ un ⎞ eλn t , (2) Затем решение системы (5) 1 λt +μ1l ⎜ ⎟ ∑ ⎝ ⎠ n =1 n ⎜ ⎟ ⎝ vn ⎠ u = aeλt +μ l , v = be 37 Математика, механика, информатика запишем в векторном виде: Зная число корней, определим сами корни, вос- пользовавшись формулой [3]: ⎛ 1 ⎞ ⎛ u ⎞ ⎜ v ⎟ = ⎜ ⎟ 1 ⎜ (λ− Aμ + D k * )⎟ μ1l +λt + 1 f ′ (λ ) ∫ λ N d λ = ∑ nk ak , ⎝ ⎠ ⎜ ⎝ 1 1 1 ⎟ ⎠ (7) 2πi f (λ ) k =1 ⎛ 1 ⎞ где ak – корень уравнения f (λ ) = 0; nk – его крат- ⎜ ⎟ 2 ⎜ (λ− Aμ + D k * )⎟ μ2l +λt . ность. ⎝ D1 2 1 ⎟ ⎠ Согласно этой формуле составим систему Для определения λ и, следовательно, λ +λ +… + λ = 1 f ′ (λ ) λ d λ, μ1 (λ ) и μ2 (λ ) воспользуемся граничными условиями 1 2 n 2πi ∫ f (λ ) (6). Подставляя решение (7) в граничные условия λ2 + λ2 +… + λ2 = 1 ∫ λ2 f ′ (λ ) d λ, при l = 0 и l = L , получим: 1 2 n 2πi C f (λ ) a eλt ⎛ λα − γ + (λβ − δ ) 1 (λ − Aμ + D k* )⎞ + . . . . . . . . . 1 ⎜ 1 1 1 1 1 1 ⎟ ⎝ 1 ⎠ λN + λN +… + λN = 1 N f ′ (λ ) λ d λ, +a eλt ⎛ λα − γ + (λβ − δ ) 1 (λ − Aμ + D k* )⎞ = 0, l = 0, 1 2 n 2πi ∫ f (λ ) 2 ⎜ 1 1 1 1 2 1 ⎟ ⎝ 1 ⎠ a eμ1L ⎛ λα − γ + (λβ − δ ) 1 (λ − Aμ + D k* )⎞ + 1 ⎜ 2 2 2 2 1 1 ⎟ ⎝ 1 ⎠ которая решается внутри достаточно большого кон- тура C. При исследовании асимптотики необходимо исхо- дить из конкретных граничных условий. Поэтому +a eμ2 L ⎛ λα − γ + (λβ − δ ) 1 (λ − Aμ + D k* )⎞ = 0, l = L. возьмем задачу, например из работы [2], и применим 2 ⎜ 2 2 2 2 2 1 ⎟ ⎝ 1 ⎠ Итак, мы имеем систему уравнений для коэффи- к ней конечно-разностную методику. Граничные условия из работы [2] приведем к виду циентов решения (7) a1 , быть равен нулю: a2 . Ее определитель должен H dv = c u − wc2 v − c v, k dt 1 V 2 l = 0, f (λ) = eμ2 L ⎛ λα − γ + (λβ −δ ) 1 (λ − Aμ + D k * )⎞ × H d c2 du = c v − L dt 2 Vd c1 u, L l = L, ⎜ 1 1 1 1 1 1 ⎟ d d ⎝ 1 ⎠ ⎛ × λα − γ + (λβ − δ ) 1 (λ − Aμ + D k * )⎞ = где H k , H d , c1 , c2 , Vd , Ld , W – известные коэф- ⎜ 2 2 2 2 2 1 ⎟ ⎝ 1 ⎠ (8) фициенты, откуда для системы (6) получим: = −eμ1L ⎛ λα − γ + (λβ −δ ) 1 (λ − Aμ + D k * )⎞ × α1 = 0, β = c1 H k γ1 = c1 , δ = −c ⎛ W ⎞ + 1⎟ , ⎜ 2 2 2 2 1 1 ⎟ ⎝ 1 ⎠ Vd H c γ = − cV ⎝ Vd ⎠ δ = c × λα − γ + (λβ − δ ) (λ − Aμ + D k * ) . α2 = d 1 , β2 = 0, 2 1 d , 2 2 . ⎜ 1 1 1 1 2 1 ⎟ ⎝ 1 ⎠ Ld Ld Вычислив корни λ этого трансцендентного урав- нения, мы тем самым найдем спектр исходной задачи, который, вообще говоря, содержит бесконечное число значений λ . Поэтому процесс решения уравнения (8) В уравнении (8) исключим младшие члены, не со- держащие λ и первого порядка относительно λ, и по- лучим уравнение для достаточно больших λ: ⎛ δ ⎞ α − 2 1 + β k разобьем на два этапа: сначала определим младшие λ ( A+ B )L ⎜ 2 D ⎜ B ⎟ ⎟ 1 корни λ внутри достаточно большого замкнутого e AB =− ⎝ 1 ⎝ ⎠ ⎠ . =β ⎛ =A ⎞ контура С и получим точное численное решение уравнения (8), а затем исследуем асимптотику при λ→ ∞ . Для получения точного решения уравнения (8) воспользуемся методами теории функций комплекс- ного переменного. Известно, что число нулей внутри α2 λ 1 ⎜1 + ⎟ D1 ⎝ ⎠ Числитель правой части, не содержащий λ , обозна- чим с. Тогда λ ( A+ B )L c e AB =− . замкнутого контура С аналитической функции f (λ ) , β1 ⎛1 =A ⎞ α2 λ ⎜ + ⎟ не имеющей полюсов, определяется по формуле [3]: D1 ⎝ B ⎠ 1 N = 2πi ∫ f ′ (λ ) f (λ ) d λ . Имея в виду, что λ – комплексное число, т. е. λ = ξ + iη, получим: 38 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева ⎛ ⎞ При исследовании полученного спектра {λn } λ ( A + B ) L ⎜ = ln ⎜ − ⎟ c ⎟ − ln λ , нужно прежде всего обратить внимание на Re λn . AB ⎜ β1 ⎛1 =A ⎞ ⎟ Если среди комплексных чисел λn имеются такие, ⎜ α2 λ ⎜ + ⎟ ⎟ ⎝ D1 ⎝ B ⎠ ⎠ для которых Re λ > 0 , то коэффициенты системы та- ⎛ ⎞ ковы, что режим будет расходящимся, и поэтому сле- ⎜ ⎟ дует изменить значения параметров. ξ= AB ln ⎜ − c ⎟ − AB ln η, 1 ( A + B ) L ⎜ β1 ⎛1 =A ⎞ ⎟ ( A + B ) L В заключение рассмотрим величину τ= max – Re λ ⎜ α2 λ ⎜ + ⎟ ⎟ ⎝ η= − D1 ⎝ 2πnAB ( A + B ) L B ⎠ ⎠ + const. время затухания самой медленной из гармоник, т. е. время установления для всей системы. От этой вели- чины зависит выбор шага интегрирования по времени Перенумеровав полученные λ в порядке возрас- тания Re λ , решение системы (3) запишем в виде суперпозиции стоячих волн: ⎛ ⎛ 1 ⎞ при расчетах по конечно-разностной методике. При расчетах периодических режимов следует иметь в виду, что существует возможность резонанса при совпадении частоты внешнего воздействия с одной из собственных частот (мнимой частью λ ). Приведем результаты расчетов, полученные с по- ⎛ u ⎞ = eλnt ⎜ a ⎜ ⎟ eμ1l + ⎜ v ⎟ ∑ ⎜ n ⎜ 1 (λ − Aμ + D k * )⎟ мощью конечно-разностной методики для прогнози- ⎝ ⎠ n ⎜ ⎜ D n 1 1 ⎟ рования поведения решения задачи (3), (4) (рис. 1–4). ⎝ ⎝ 1 ⎠ ⎛ 1 ⎞ ⎞ Были использованы следующие коэффициенты: B ⎜ ⎟ n ⎜ (λ − Aμ + D k * )⎟ μ2l ⎟ ⎟ A = 36 м/с, B = 360 м/с, D1 = D2 = k y V , ⎝ D1 2 1 ⎟ ⎟ ⎠ ⎠ V = 73, 07 кмоль/ч, k * = V ⋅ A , L ⋅ B α= 45, 21 кмоль/ч, где an , bn выбираются таким образом, чтобы удов- H k = 30 кмоль/ч, H d = 50 кмоль/ч, летворить заданному начальному условию. W = 76, 59 кмоль/ч. u 1 0,9 0,8 0,7 0,6 0,5 0,4 0,3 0,2 u(40,t) u(20,t) u(0,t) 2 4 6 8 10 12 14 16 18 20 22 t, ч 1 2 3 0,1 u(0,t) u(20,t) u(40,t) Рис. 1. Кривые переходных процессов при возмущении по концентрации вверху колонны для точек l = 0, 20, 40 м: 1′, 2′, 3′ – при увеличении u ( L,0 ) на 50 %; 1, 2, 3 – при уменьшении u ( L,0 ) на 50 % 39 v 0,7 0,6 0,5 Математика, механика, информатика u(40,t) u(20,t) u(0,t) 2′ 1′ 3′ u(40,t) 3 2 4 6 8 10 12 14 16 18 t, ч 0,4 0,3 u(20,t) 1 u(0,t) Рис. 2. Кривые переходных процессов при возмущении по концентрации вверху колонны для точек l = 0, 20, 40 м: 1′, 2′, 3′ – при увеличении v ( L,0) на 50 %; 1, 2, 3 – при уменьшении v ( L,0) на 50 % u 1 0,9 0,8 0,7 0,6 0,5 t = 2 t = 6 t = 10 t = 14 t = 20 10 20 30 40 l, м Рис. 3. Профиль концентраций в жидкости по длине колонны при возмущении на 50 % вверху колонны для различных t, ч v 0,6 0,5 0,4 t = 2 t = 6 t = 10 t = 20 10 20 30 40 l, м Рис. 4. Профиль концентраций в паре по длине колонны в различные моменты времени при возмущении на 50 % вверху колонны для различных t, ч 40 Вестник Сибирского государственного аэрокосмического университета имени академика М. Ф. Решетнева Для этих коэффициентов все Re λ< 0 , что обеспе- чивает сходимость процесса. Кроме того, можно оп- (ректификационных установок), основанный на мате- матическом аппарате, содержащем дифференциаль- ределить время установления τ= max 1 Re λ = 6, 39 ч, ные уравнения в частных производных. Для исследо- вания процессов установления применялся метод стоящих волн, который позволяет получить сведения за которое амплитуда самой медленной гармоники уменьшается в с раз. В сложных химико-технологических установках, состоящих из ряда соединенных аппаратов и имею- щих систему потоков взаимодействия их сред, пара- метры которых чувствительны к возмущениям, все более важным становится построение высокоэффек- тивных систем контроля и управления. Однако преж- де чем создавать такие системы управления, полагая основным ее звеном управляемый объект, необходи- мо улучшить статические и динамические характери- стики этого объекта. В данной статье предложен подход к решению задачи моделирования нестационарных режимов объ- ектов с рециркуляцией взаимодействующих потоков не только о качестве и характере процессов установ- ления, но и определить важные характеристики вре- мени установления, резонансные свойства системы и выявить области неустойчивости в пространстве па- раметров системы