CC BY
23
Серия «Математика»
2012. Т. 5, № 3. С. 2-17
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 517.97
Задача разделения смесей в ректификационных колоннах *
М. А. Аргучинцева
Иркутский государственный университет
В. П. Поплевко
Иркутский государственный университет
Аннотация. В данной статье исследована прикладная задача управления ректификационной колонной. Процесс предназначен для получения из многокомпонентной смеси продуктов с заданной концентрацией. Математическая модель процесса ректификации представляет собой систему гиперболических уравнений первого порядка. Граничные условия системы определяются из систем обыкновенных дифференциальных уравнений.
Ключевые слова: гиперболические уравнения, численные методы, ректификационная колонна
Интересным приложением задач оптимального управления является моделирование химических и технологических процессов в промышленности. Например, процесс ректификации широко применяется в нефтегазоперерабатывающей, химической, нефтехимической и других отраслях. Он предназначен для получения из многокомпонентной смеси продуктов с заданной концентрацией. Основу процесса ректификации составляет теплообмен, массообмен и гидродинамика взаимодействующих потоков. Этот процесс характеризуется большим числом параметров, связанных между собой сложными зависимостями, и описывается системами дифференциальных уравнений в частных производных [1]. Значительная часть параметров является функциями временной и
Введение
* Работа выполнена при финансовой поддержке ФЦП «Научные и научнопедагогические кадры инновационной России на 2009-2013 годы» и РФФИ (грант 11-01-00713).
пространственной координат. Цель управления процессом ректификации - поддержание заданного состава целевого продукта. Поставленная задача принадлежит классу задач оптимизации гиперболических систем с управляемыми дифференциальными связями на границе, где управления выбираются из класса гладких функций и удовлетворяют поточечным или интегральным ограничениям.
Поскольку на каждой итерации оптимизационного алгоритма приходится интегрировать сложную начально-краевую задачу для системы гиперболических уравнений, поэтому в данной работе авторы более подробно остановились на методе ее исследования.
1. Постановка задачи
Рассмотрим математическую модель процесса разделения смесей в ректификационных колоннах [1]:
д(НХ д(ЬхЛ = - *) + *«,
dt ds
д(НУ y) + d(M_ k (y*_ у.) + ф
dt + ds = kyi(Vi Vi) + ФУ1,
N N
Exi = ^Vi = 1, i = 1,N.
i=1 i=l
Здесь xi = xi(s,t), Vi = Vi(s,t) - концентрации i-го компонента в жидкой и паровой фазах; L = L(s,t), V = V(s,t) - потоки жидкости и пара в колонне; Hx = Hx(s,t), НУ = Hy(s,t) - удерживающие способности колонны по жидкости и пару; ФХ1 = ФХ1 (s,t), ФУ1 = ФУ1 (s,t) -плотности вводимых потоков i-го компонента исходной смеси в жидкой и паровой фазах; s € [so,si] - координата вдоль ректификационной
колонны; t € [to,ti] - время. Коэффициент массопередачи в паровой
фазе имеет вид:
кУг = кУ = kV(s, t), k = const.
Концентрация паровой компоненты в равновесной ситуации находится по формуле:
V*(s,t) = p(s,t)xi(s,t),
где p)(s, t) определяется из эксперимента.
Вверху и внизу колонны имеются две емкости - конденсатор (”d”) и испаритель (”k”), в которых накапливаются соответствующие "легкие" и "тяжелые" компоненты разделения исходной смеси. Исследуем подробнее процессы массообмена, проходящие в испарителе. Поток
ректификационной жидкости Ь^(Ь) = Ь(во,Ь), выходящий из нижней части колонны, поступает в испаритель (в = во). Часть потока испаряется в кипятильнике и возвращается в колонну в виде пара Vк(і) =
V(во, і). Другая часть отбирается в виде готового продукта Ш(і). Таким образом, концентрация компонентов жидкости Хкі(і), находящейся в испарителе, определяется из уравнения материального баланса [1]
^Нхк= Ък(і)хі(во, і) - Vk(І)Уі(во, і) - Ш(і)хкі(і),
Хкі(іо) = Хіо(во), і = 1,Ы.
Здесь Уі(во,і) находится из дополнительного соотношения, учитывающего различные типы кипятильников, стоящих в испарителе,
Уі(во,і) = (уі (во,і) - Хкі(і))ак + Хкі(і),
ак = 0 - для полного испарителя и ак = 1 - для парциального испарителя. Общее количество жидкости в испарителе Нхк(і) рассчитывается из уравнения
^к(і) - ^(і) - Ш (і), Нхк(іо) = Нхко.
Рассмотрим уравнения массообмена в конденсаторе (в = ві). Паровой поток Ці(і) = V(ві,і), выходящий вверху колонны, поступает в конденсатор. Часть потока жидкости, выходящего из конденсатора, отбирается в виде готового продукта О(і), другая часть поступает в верхнюю часть колонны в виде потока орошения Ь^(і) = Ь(ві,і). Уравнение покомпонентного материального баланса для концентрации жидкости в конденсаторе х^(Ь) [1]:
= Vd(і)ydi(і) - (Ьа(і) + Б(і))хаі(і),
Хсіі(іо) = Хіо(ві), і = 1,М,
где
Уdi(і) = Уі(ві,і)+ ad(у*(ві,і) - уі(ві,і)),
ad - эффективность конденсатора (ad = 0 - для полного конденсатора и ad = 1 - для парциального конденсатора). Количество жидкости в конденсаторе Нх^(і) определяется из уравнения общего материального баланса потоков:
= Vd(і) - Ld(і) - Б(і), Нх^(іо) = НхМ-аі
Введем основные предположения:
1) подвод сырья осуществляется только в жидкой фазе в ст-окрестности точки s = s*
ф№ (s,t) = о,
$xi(s,t) = XftГх(г)фх1 (s), s* - a < s < s* + a,
где Xfi - концентрация компонента в жидкости потоков сырья, Fx(t)
- вводимый поток сырья, фх^ (s) - функция распределения потока по длинне колонны - заданные функции;
2) в ректификационной колонне с увеличением потоков L(s, t) и
V(s,t) увеличиваются удерживающие способности Hx(s,t) и Hy(s,t). Будем считать эту зависимость прямо пропорциональной
L(s,t) = С1 = const, ^^7^ = С2 = const;
Нх(в,і) ’ Ну (в,і)
3) поток пара в колонне зависит только от времени і и равен потокам пара в конденсаторе и испарителе:
V(в, і) = V(і), V(і) = Vd(і) = ^к(і);
4) потоки жидкости в колонне, конденсаторе и испарителе зависят только от времени і и определяются следующим образом:
L(в, і) = L(і) + L*(і),
^х(і), во < в < в* О,
З*+(7
L*(і) = ^ Рх(і) / Фхі(в) (їв, в* - О < в < в* + О,
З* —О
0, в* + о < в < ві,
Ld(і) = L(і), Lk (і) = L(і) + ¥х(і);
5) должно выполняться следующее условие между количеством исходной смеси, вводимой в колонну, и потоками готовых продуктов, отбираемых в конденсаторе и испарителе: Ех (і) > О (і) + Ш (і). Если колонна работает в "статическом" режиме (количество исходной смеси равно количеству готового продукта), справедливо строгое равенство в каждый момент времени і: Рх(і) = О (і) + Ш (і) или за весь период работы колонны:
*1
J[Fx(t) _ D(t) _ W(t)] dt = О;
to
6) рассматривается случай полного испарителя (слк = 0) и конденсатора (ad = 0), тогда
Хкі(і') уі (во,і'),
ум(г) = Уг(вг,г), хЛг(г) = Хг(в1,г).
С учетом всех предположений система уравнений в частных производных примет вид:
дХгд^’ ^ -С19Хгд1 ’ ^ = А1(8,г)хг(8,г) + В1(8,г)уг(в,г) + Р1(в,г), ду%((г ^ + С2 9Уг(д88’г) = А2(8’ г)хг(8, г) + В2(8, г)уг(8, г) + Г2(8’ г), (1)
где
л ( ^ _ С1Гх(г)фхг(8) - слкУ(г)р(8,г) - (ь(г) + ь*(г))'
А1(8’г) =---------------
ь(г) + ь*(г)
А2(8’г) = С2кр(8’г)’ Г1(8’г) = Сц£)(+фщ8’ Ы8’г) = 0’ п / ч С1кУ(г) у'(г).
В1 (8-‘) = тщтпщ’ В'2(8Л) = -(кС2 + т)'
N N
Х
г=1 г=1
ЕХг(8’г) = 1’ ^Уг(8’г) = 1’ Хг(8’ г )> 0’ Уг(8’г)> 0’
80 < 8 < 81 ’ г0 < г < г1’ г = 1’Ы.
Начальные условия имеют вид:
Хг(8’го )= Хг0 (8)’ Уг(8’го) = Уг0 (8)’ 80 < 8 < 81’ г = 1’Ы. (2)
Граничные условия при го < г < г1 на границе 8 = 8о:
й(нхк((уг(8о’)) = [Ь(г) + Рх(г)-]Хг(80’г) - [у(г) + ш(г)]Уг(80’г)’
Уг(80’ го') — Уг0 (80) ’ г — 1’ N; (Шхк (г)
(г
= ь(г) + Гх(г) - у (г) - ш (г)’ НхМ = Нхк0; (3)
и на границе 8 = 81:
((Иха(г)Хг (81’г)) (г
= у(г)Уг(81’г) - [ь(г) + о(г)]Хг(81’г)’
Хг(81’г0) = Хг0(81)’ г = 1’М;
(Нхм(г)
(г
= у (г) - ь(г) - Б(г)’ НхА (г0) = Нхю. (4)
Параметры задачи с1, с2, к, 80, 81, 8*, г0, г1, а, Нхк0, Нх^0, р>(8’г), Fx(г), фх1 (8), Х/г, Хг0(8), Уг0(8) считаются заданными.
Поставим задачу оптимального управления процессами разделения смеси в ректификационной колонне, описываемую гиперболической системой первого порядка (1) с начальными условиями (2) и граничными условиями (3), (4). В качестве управлений здесь выбираются потоки готовых продуктов, отбираемых в испарителе Ш(г) и конденсаторе О(г). Особенностью является то, что управления входят в правые части дифференциальных связей (3), (4) на границах потока 8 = 80 и 8 = 81. Предполагается, что Ш(г) и Б(г) принадлежат классу гладких функций.
В качестве целевого функционала 3 задается суммарное отклонение концентраций в выходных потоках Хмг(г)’Хкг(г) от заданных значений
вц’ $2г’ г = 1’М:
N *
3 = Е / [Ки(Хаг(г) - Ои)2 + К2г(Хкг(г) - 02г)2} йг ^ тгП’
г=1*0
который в случае полного испарителя и конденсатора принимает вид
N *
3 = Е / [К1г(Хг(81 ’г) - Он)2 + К2г(Уг(80’ г) - в2г)2] (г ^ ШЪП. (5)
г=ч0
Здесь Кц’ К2г - весовые коэффициенты, определяющие ценность продукта.
Таким образом, задача (1)-(5) принадлежит классу задач оптимизации гиперболических систем с управляемыми дифференциальными связями на границе, где управления выбираются из класса гладких функций и удовлетворяют поточечным (интегральным) ограничениям. В работе [1] данная задача рассматривалась в классе кусочно-непрерывных функций. Методика [2,3] позволяет исследовать ее в классе гладких управляющих функций.
На каждой итерации оптимизационного алгоритма приходится интегрировать начально-краевую задачу для системы гиперболических уравнений (1)-(4), поэтому более подробно остановимся на методе ее исследования.
2. Метод исследования
Одним из наиболее эффективных методов исследования систем гиперболических уравнений является метод характеристик [4]. Он заключается в отыскании кривых, называемых характеристиками, вдоль ко-
торых уравнения в частных производных превращаются в обыкновенные дифференциальные характеристические уравнения. Отличительной особенностью метода, по сравнению с другими разностными алгоритмами, является минимальное использование операторов интерполирования и, в связи с этим, максимальная близость области зависимости разностной схемы и области зависимости системы дифференциальных уравнений.
В области Б х Т, (Б = ], Т = [Ь0,Ь1]) построим характе-
ристическую разностную сетку, узлы которой являются пересечениями характеристик первого (8(1)) и второго (8(2)) семейств, задаваемых уравнениями:
ёв(1) ёв(1)
М Сі, М °2'
В общем случае характеристическая сетка является криволинейной, что существенно затрудняет численную реализацию метода характеристик из-за сложной логики построения фронта расчета. Однако в данном случае из физической специфики задачи о ректификационной колонне (1)—(4) вытекает, что производственные коэффициенты
С1 = СОП8Ь, С2 = СОП8Ь.
Следовательно, характеристики обоих семейств представляют собой прямые линии, и можно построить разностную сетку, свободную от указанного недостатка.
Разобьем отрезок Б = [80,81] на т равных частей [8у,8^+1], положив 8 у = 80 + і А8, і = 0,1,...,т; А8 = 31—30. Через точку 8 = 81 проведем нулевую характеристику второго семейства 8(2). На этой характеристике отметим точки Му (8у ,щ), где координаты у вычисляются по правилу:
8і — 81
щ = Ьо + —-------, І =0,1,...,т
С2
Из точек Му, і = 0,...,т выпускаем характеристики первого семейства 8(1) и найдем точки их пересечения с осью Ь = Ьо'
Я1з = 8у + СлЩ — Ьо), і = 0,...,т,
8о < Я1У < 81.
В точках q1j, і = 0,...,т задаются начальные условия (2) для характеристик обоих семейств.
Из точек пересечения характеристики первого семейства 8(1) с прямой 8 = 8о'
Ьі0 = ио + іАЬ, АЬ = (— + —)А8, і = 0,1,...,Е
С1 С2
проводим характеристики второго семейства 8(2). Здесь К - номер последней характеристики второго семейства, попадающей в допустимую область Б х Т. Точки пересечения характеристик 8(2) с прямой 8 = 81 имеют вид:
Ьіт — щт + iАt, і — ° 1 ,...,П.
В точках пересечения характеристик второго семейства с прямой Ь = Ьо'
^2ї = 81 + С2(Ь0 — Ьіт), і = 0 1,...,К,
80 < q2i < 81
задаем начальные условия (2).
Точки пересечения характеристик первого и второго семейства определяются по формулам:
Ьч = Щ + іАЬ, і = 0,1,...,К, і = 0,1,...,т, є Б х Т.
Таким образом получили характеристическую разностную сетку с узлами (8у, Ьіу), і = 0,1,...,К, і = 0,1,...,т.
Введем следующие обозначения:
1) значение функции в узле характеристической сетки 8,Ьіу) будем обозначать индексами [.. .](8у, Ьіу) ~ [.. .]і,ч ,
2) р1 = ЬіЧ Ьі—1,Ч+1 , р2 = Ьіу Ьі,Ч — 1 ,
3) Аь = Ьіу Ьі—1,ч.
Выведем расчетные формулы для внутренних точек области Б х Т. Первое из уравнений (1) интегрируется по формуле трапеций вдоль характеристики первого семейства (8(1)) в направлении Ь—1ч+1 ^ , а
второе - вдоль характеристики второго семейства (8(2)) в направлении Ьіч—1 ^ Ьіу. Тогда получим следующую неявную разностную схему:
Хі,Ч Хі—1,Ч + 1 1 ( Л І О I Е1 I
-------------- 7л {А1ічхі,Ч + В1і,ЧУі,Ч + Рич +
р1 2
+А1і—1,Ч+1 хі—1,ч+1 + В1 і—1ч+1уі—1ч+1 + Fli—l^j+l},
Уч----~ = о {А2і,ухі,Ч + В2ічуі,Ч + Р2і ч +
Р2 2
+А2іЧ—1хіч — 1 + В2іч — 1уіч — 1 + Р2іч — 1}.
Введем обозначения:
Щ+\ = хі—1,Ч+1 + р~ {А1і—1,Ч+1 хі—1,Ч+1 + В1 і—1ч+1уі—1ч+1 + Fli—l,j+l},
Р2
Я)—1 = уі,У—1 + {А2іЧ—1хі,У—1 + В2іч — 1уіч — 1 + Р2іч — 1}.
Тогда наша система примет вид:
хі
Р1
Р1
2
2
2
У* ,3 = Р^ А2* >3 Х* >3 + Р2 В2* >3 У* >3 + ( ^2 ^ ^
Рассматривая данные выражения для неявной разностной схемы как систему двух линейных алгебраических уравнений относительно х*,3 и У
і, ч-
(1 — у А1і ч )хі ч — Р^ В1і ч уі ч = У Р1і ч + Бч+1, —Р2 А2і ч хі ч + (1 — Р22 В2і ч )уі ч = Р22 ^ ч + Я)—1,
2
2
можно получить явные формулы перерасчета разностного решения на г-ом слое по формулам Крамера:
Ах
Д.
і ч
А
уі ч
уг]
А
где
А
(1 — Р1 Аич) — Рт Вцч
—РТ А2 ч (1 — Р2 В2і,ч)
гз
Ах
(Рт р1і ч + Л)
( РТ Р2і ч + Яч—1)
-Рт В1 ■
2 1 і
і ч
(1 — РТ В2 іч )
г3
А
Угз
(1 — РТ Аіі ч) (РТ Рі і ч + Б—) — РТ А2іч ( РТ Р2іч + Я) — 1)
Окончательно с учетом обозначений:
М3 = {Ви3 Р2*3 - В2’3 РЧ3 ) + Т Р1 *3 ],
мч = Кіч [—Р1р2 (Аіі ч Р2і,ч — А2іч Рііч ) + Ц. Р2іч ],
= Кіч [(1 — ^В2і,ч) VР1і— 1 ч+1 + VВііЛР2іч — 1]
Р2
Р1
р1 о р2
N2 = Кіч [^ А2іч 2 Р1і— 1 ч+1 + (1 — 2 Аіі ч )'2 Р2іч — 1], БіЧ = Кч(1 — РТ В2іч)(1 + Рт А1 і— 1 ч+1),
г Р2 л Р1
Р1
, Р2
Би = Кч^Т А2 іч(1 + V А1 і— 1 ч+1),
Р2
Р1
ЗАДАЧА РАЗДЕЛЕНИЯ СМЕСЕЙ ІІ
Vj = Rij (1 _ f B2ij if Bn- j + L
V( 2) ___ R P2 A Pl B
Vij = Rij Y A2iJ ~2 Bu-1 j+1,
Til = Rij 2 Blij f A^j-L Tj = Rij (1 _ У Alij if A^j-l, wj = Rij f Blij (1 + P2 Buj-l), w(2 = Rij(1 _ f Alijji(1 + ^2B2j-li
получим явные формулы перерасчета разностного решения на i-м слое по заданным значениям на (i _ 1)-м слое:
Xi,j = Si,j Xi-1 j+l + Ti,j Xij-1 + Vi, jlyi-1 j+l + Wi,j yij-1 + M(j + Ni,j ,
Уі j = S%Xj-l j+l + Tj, jjxj j-1 + Vj У j-1 j+l + Wjjjyi j-1+Mj +Nj (6)
для внутренних точек (sj, tij) области S x T. Предлагаемая разностная схема (6) является устойчивой и обеспечивает второй порядок аппроксимации.
Получим расчетные формулы на границе s = s, , где должны выполняться условия (4). Рассмотрим два случая:
І. Удерживающая способность конденсатора Hxd = Hxd0 = const, т.е. рассматривается режим работы ректификационной колонне, при котором количество жидкости в конденсаторе постоянно. Тогда условия (4) перепишутся в виде:
dx(sl,t) _ L(t) + D(t\(y(sl,t) _ x(sl,t)), x(sl,to) = xo(sl).
dt H
хй0
2. Удерживающая способность конденсатора меняется с течением времени Нха = Нх^(Ь). Тогда
йх(81,Ь) V(Ь) , , л , ,, , . , .
Л------ = ГТ (Лу(81,Ь) — х(81,Ь)), х(81,Ьо) = хо(81),
М Нхй(Ь)
Н^ = V (Ь) — Ц(Ь) — Б(Ь), Нха(Ьо) = НхМ.
К границе 8 = 81 подходят характеристики второго семейства 8(2), на которых:
йу(81,Ь)
dt
= A2(Sl,t)x(Sl,t) + B2(sl,t)y(s,t) + F2(Sl,t).
Запишем, используя формулу трапеций, разностные аналоги этих уравнений в случае Hxd = Hxd0 = const:
xi,m xi—1,m — {Kdi,m(yi,m xi,m) + Kdi—1,m(yi—1,m xi—1,m)},
2
yi,m yi,m—1 = f {A2i,m—1Xi,m—1 + B2i,m—1yi,m—1 + F2i,m—1 + +A2i,mXi,m + B2i,myi,m + F2i,m}, i 0, 1,...,n,
где
Li,m + Di,m Li—1,m + Di—1,m
Kdi,m 77 , Kdi—1,m 77 ,
Hxdo Hxdo
Перепишем последние выражения в более удобном виде:
2 2 2
(_ Y A2 i,m )xi,m + (1----2 B2i,m')yi,m = (f A2i,m-1)xi,m-1 +
2 2 + (1 + Y B2i,m—1)yi,m—1 + ~Y (F2i,m— 1 + F2i,m),
, At . , At . , At
(1 + 2 Kdi,m)xi,m + ( 2 Kdi,m)yi,m = (1 ^ Kdi—1,m)xi—1,m +
AtK
+ 2 Kdi—1,myi—1,m.
Решим систему линейных алгебраических уравнений относительно Xj,m и yi}m. Тогда с учетом обозначений:
Gd^l1(i) = Ad(i)( Pf B2 i,m _ 1)[1---Y~Kdi-1,m],
Gd^2^(i) = Ad(i) P2 A2i,m(1-----2 Kdi—1,m),
PdW(i) = Ad(i)( 2 B2i,m _ 1) AfKdi-1,m,
Pd(2)(i) = _Ad(i)P2AKdi-
^d( ^(i) — Ad(i) a Kdj,m(F2j,m—1 + F2i,m^
Ud^(i) = Ad(i) -2- (F2j,m + F2j,m-l)(1 + AГKdi,m),
Ш(1(і) = Ad(i) Y~Kdi,m(1 + — B2i,m— 1),
= Ad(i)(1 + P^ B2i,m—1)(1 + YГKdi,m),
№d( ^(i Ad(i) 4 A2i,m—1Kdi,m
№d(2 (i) = _Ad(i) P2 A2i,m—1(1 + A^Kdi,m)
получили явные расчетные формулы на границе s = s,:
Xj,m=№d('l) (i)Xj,m-1+nd('l) (i)yi,m-1+Gdi'l) (i)Xj-l,m + Pd!'1 (i)y—1,m + Ud^1 (І), yi,m=^d{2) (i)Xj,m-1+nd{'2) (i)yi,m-1+Gd('l) (i)Xj-l,m + Pd^ (i)y—1,m + Ud^ (i).
(7)
Здесь x0,m, y0,m известны из начальных условий (2).
Аналогичная процедура проводится и в случае Hxd = Hxd(t), для которого будут справедливы формулы (7) с точностью до значений коэффициентов:
Vi,m Vi—1,m
Kdi,m H , Kdi—1,m H ,
Hxdi Hxdi—1
Hxdi — Hxdi—1 + 2 {Vi,m + Vi—1,m Li,m Li—1,m Di,m Di—1,m}.
Рассмотрим условия в испарителе (s = so). Возможны два случая:
1. Удерживающая способность испарителя Hxk = Hxk0 = const. Тогда условия (З) примут вид:
dy(d°,t') = V (t)H W (t) (x(so,t) _ y(so,t)), y (so ,to) = yo(so). at Hxko
2. Удерживающая способность испарителя меняется с течением времени Hxk = Hxk (t). Тогда
dy(so ,t) L(t) + Fx(t) ,,,,,,, , .4 f s
—у,— =---------Й-(x(so,t) _ y(so,t)), y(so,to) = yo(so);
dt Hxk(t)
= L(t) + Fx(t) _ V (t) _ W (t), Hxk (to) = Hxko.
К границе s = So подходят характеристики первого семейства s(1), на которых:
dx(s0,t)
dt
= Al(so,t)x(so ,t) + Bl(so,t)y(so,t) + Fl(so,t).
Запишем разностные аналоги этих уравнений, используя формулу трапеций:
уг,0 — уг-1,0 = {Ккг,0 (хг,0 — уг,0) + Кк—1,0 (хг-1,0 — уг-1,о)},
Хг,0 — Хг-1,1 = р- {А1 г,0Хг,0 + Е1г,0уг,0 + Р1г,0 + +А1г-1,1Хг-1,1 + Е1г-1,1уг-1,1 + Рц-1,1},
где
^ _ Уг,0 + ^г,0 ^ _ Уг-1,0 + №г-1,0
Ккг,0 = ----77-) Ккг-1,0 = ----------------------77-•
нхк0 нхк0
Перепишем в более удобном виде:
(1 — ^2 А1г,0 )хг,0 + ( — р В1г,0)уг,0 = (1 + р" А1г-1,1)Хг-1,1 + + р" Е1г-1,1уг-1,1 + р (Р1 г,0 + P1г-1,l),
АЬ А1 А1
---2 Ккг,0Хг,0 + (1 + у Ккг,0)уг,0 = ~ Ккг-1,0(хг-1,0 +
+(1-----2'Ккг-1,0)уг-1,0 •
С учетом обозначений:
Ок {1)(г)= вг Р~А~Ец,0Ккг-1,0,
Ок(2)(г)= вг А (1 - ^ А1г,0)Ккг
2 2 — кг,0!
Рк {1)(г)= вг Р2 Еи>0(1 -А±Ккг-1,0),
Рк(г) = вг(1 — А1г,0)(1------АКкг-1,0))
Шк(1)(г) = вг(1 + А Ккг,0) ^2 (Р1г,0 + P1г-1,l),
Шк(22(г) = вг — Ккг,0(Р1г,0 + Р1г-1,1))
Пк {1)(г)= вг(1 + А Ккг,0) 2 Е1 г-1,1,
Пк[2)(г)= вг ^ В—
4 Е1г-1,1Ккг,0,
Ик {1)(г) = вг[(1 + А Ккг,0)(1 + 2 Ан-1,1)],
№ {2)(г)= вг А КЧ0(1 + р2 Ан-1,1)• расчетные формулы имеют следующий вид:
Хг,0 = Ок {1)(г)х г-1,0 +Рк {1](г)у г-1,0 + Пк ^(г)уг-1,1 +^к(1'}(г)хг-1,1 +Шк <'1\г),
;,о = Gk ^ (i)x—i,0+Pk ^ (i)y—i,o+Пк ^ (i)y-i,i +^к ^ (i)xi-\,i +ши ^ (i),
i = n*,...,R, (8)
где п* - номер первой характеристики второго семейства, пересекшей прямую в = в0.
Аналогично рассматривается случай Нхк = Нхк(Т) для системы разностных уравнений, где
^ _ Рг,0 + Рхг,0 ^ _ рг-1,0 + Рхг-1,0
Ккг,0 = ------------Н-! Ккг-1,0 = -Н----------•
Нхкг Нхкг-1
Нхкг=Нхкг-1 + А^Т {Рхг,0+Рхг-1,0+рг,0+рг-1,0 — Уг,0 — Уг-1,0 — '^г,0 — '^г-1,0}•
Таким образом, и в этом случае расчетные формулы на границе в = в0 имеют вид (8).
Для проверки правильности построения разностной сетки и эффективности введенных конечно-разностных схем, решим несколько тестовых задач, сравнивая полученные решения с аналитическими.
Пример 1.
Рассмотрим задачу:
^ - «1 ^ = ^ у(в',)'
0^ + «2^ хМ + _+_у(аЛ
«1 = 1,«2 = 3
с граничными условиями:
йу(в0,Т) 2совТ
dt cost — 2sint dx(sl,t) —esi sint
— (x(so,t) — y(so,t)),
(y(si,t) — x(si,t))
dt s + 2)sint — esi cost и начальными условиями:
x(s, 0) = es, y(s, 0) = 0.
Аналитическое решение поставленной задачи имеет вид:
x*(s, t) = escost, y*(s, t) = (s + 2)sint.
Предложенный характеристический численный алгоритм реализован в системе MatLab 7.0. Сравнение результатов расчетов x(s, t), y(s, t)
y
с аналитическим решением х*(в,Ь), у*(в,Ь) показало, что даже на достаточно крупной сетке при в € [0, 2], Т е [0, 4], т = 60 погрешности численного решения имеют вид: тах \х* — х\ = 0, 00056; тах\у* — у\ =
0, 00027-
Пример 2.
Дана задача:
дх(в,Т) дфЛ = дТ - С1 дв = -С1х(в,Т) - в + 2У(в,t),
ду(в,г) ду(в,Т) в + 2 «2 . .
~дТ + «2^^ = — Х{в'Т) + ~2 ^
С1 = 1, «2 = 3^
Граничные условия:
d(Нхк(Т)у(в0, Т))
dt dHxk (t)
= 2costx(s0,t) _ (4cost + sint)y(s0,t), = _2cost _ sint, Hxk (0) = 1;
d(Hxd(t)xd(t)) _ sl
dt
dHxd(t)
dt
Начальные условия:
dt
esl sinty(s, ,t) _ (2esl sint + (s, + 2)cost)x(s1, t),
= _eslsint _ (s, + 2)cost, Hxd(0) = esl
x(s, 0) = es, y(s, 0) = 0.
Аналитическое решение задачи:
x*(s, t) = escost, y*(s, t) = (s + 2)sint.
Результаты расчетов при s £ [0, 2], t £ [0, 4], m = 60 дали следующие погрешности численного решения: max\x* —x\ = 0, 00063; max\y* —y\ =
0,00031.
Список литературы
1. Демиденко Н. Д. Моделирование и оптимизация систем с распределенными параметрами / Н. Д. Демиденко, В. И. Потапов, Ю. И. Шокин. - Новосибирск : Наука, 2006. - 551 с.
2. Аргучинцев А. В. Оптимальное управление гиперболическими системами / А. В. Аргучинцев. - М. : Физматлит, 2007. - 168 с.
3. Аргучинцев А. В. Оптимальное управление граничными условиями гиперболической системы на примере задачи химической ректификации / А. В. Аргучинцев, В. П. Поплевко // Тр. XV Байк. междунар. школы-семинара «Методы оптимизации и их приложения», п. Листвянка, оз. Байкал, 23-29 июня 2011 г. - Иркутск : РИО ИДСТУ СО РАН, 2011. - Т. 3. - С. 36-40.
4. Рождественский Б. Л. Системы квазилинейных уравнений и их приложения к газовой динамике / Б. Л. Рождественский, Н. Н. Яненко. - М. : Наука, 1978. -686 с.
M.A. Arguchintceva, V.P. Poplevko
A problem of fractionization process in a tower
Abstract. In this paper a process of fractionization in a tower is considered. This process is described by a system of first-order partial differential equations. The numerical experiment is carried out.
Keywords: hyperbolic partial differential equations, computational methods, fractioni-zation in a tower
Аргучинцева Маргарита Александровна, кандидат физико-математических наук, доцент, Институт математики, экономики и информатики, Иркутский государственный университет, 664003, Иркутск, ул. К. Маркса, 1 (marguch@math.isu.ru)
Поплевко Василиса Павловна, кандидат физико-математических наук, доцент, Институт математики, экономики и информатики, Иркутский государственный университет, 664003, Иркутск, ул. К. Маркса, 1 (vasilisa@math.isu.ru)
Arguchintceva Margarita, associate professor, Irkutsk State University,
1, K. Marks St., Irkutsk, 664003 (marguch@math.isu.ru)
Poplevko Vasilisa, associate professor, Irkutsk State University, 1, K. Marks St., Irkutsk, 664003 (vasilisa@math.isu.ru)
