ИССЛЕДОВАНИЕ СТАТИЧЕСКИХ И ДИНАМИЧЕСКИХ РЕЖИМОВ ПРОЦЕССОВ МАССООБМЕНА

Полный текст

Для математического описания процессов массо- обмена в ректификационных колоннах широко ис- пользуются системы уравнений в частных производ- ных [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 ч, ные уравнения в частных производных. Для исследо- вания процессов установления применялся метод стоящих волн, который позволяет получить сведения за которое амплитуда самой медленной гармоники уменьшается в с раз. В сложных химико-технологических установках, состоящих из ряда соединенных аппаратов и имею- щих систему потоков взаимодействия их сред, пара- метры которых чувствительны к возмущениям, все более важным становится построение высокоэффек- тивных систем контроля и управления. Однако преж- де чем создавать такие системы управления, полагая основным ее звеном управляемый объект, необходи- мо улучшить статические и динамические характери- стики этого объекта. В данной статье предложен подход к решению задачи моделирования нестационарных режимов объ- ектов с рециркуляцией взаимодействующих потоков не только о качестве и характере процессов установ- ления, но и определить важные характеристики вре- мени установления, резонансные свойства системы и выявить области неустойчивости в пространстве па- раметров системы

Об авторах

Николай Данилович Демиденко

СКТБ «Наука» КНЦ СО РАН

Email: secretary@sktb.krsn.ru.
– доктор технических наук, ведущий научный сотрудник СКТБ «Наука» КНЦ СО РАН. Окончил Томский государственный университет в 1963 г. Область научных интересов – моделирование, управление и оптимизация систем с распределенными параметрами.

Список литературы

Дополнительные файлы