CC BY
44
УДК 517.518
О ДЕКОМПОЗИЦИИ РАЗНОСТНЫХ СХЕМ ПРИ ЧИСЛЕННОМ РЕШЕНИИ
ДИФФЕРЕНЦИАЛЬНО-АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
В.Ф. Чистяков, Э.А. Таиров, Е.В. Чистякова, А.А. Левин
Рассматриваются квазилинейные системы обыкновенных дифференциальных уравнений (ОДУ), с тождественно вырожденной матрицей перед производной искомой вектор-функции и разностные схемы, применяемые для их решения. В работе обсуждаются условия, обеспечивающие на каждом шаге вычислительного процесса возможность последовательного решения алгебраических (конечных) уравнений и подстановки этих решений в уравнения динамики. Приведены результаты численных экспериментов для систем ОДУ, описывающих прямоточную котельную установку.
Ключевые слова: дифференциально-алгебраические уравнения, индекс, декомпозиция, разностные схемы, математические модели, прямоточные паровые котлы.
Введение
Рассмотрим вырожденную систему обыкновенных дифференциальных уравнений
A(x, t, v)x + B(x, t, v) = 0, t € T = [a, ß], (1)
где A(v, t, v) — (пхп)-матрица, B(v, t, v) —n-мерная вектор-функция, определенные в области V = Rn хТxN, x = x(t, v)—искомая вектор-функция, v—векторный параметр из множества N = [—vi, vi] х [—vi, vi] x ■ ■ ■ x [—vm, vm] c Rm, x = dx/dt.
Предполагается, что характер вырождения определяется соотношениями
det A(v, t, v) = 0 V (t, v, v) € V, (2)
и заданы начальные данные
x(a, v) = a(v), (3)
где a(v) - заданная вектор-функция, v € N.
Под решением системы (1) на отрезке Т£ = [а, а + е] С Т при значении параметра v = v* мы будем понимать любую вектор-функцию x(t,v*) € C1(T£), которые обращают исходные уравнения в тождество на Т£.
В частности, в виде системы (1), удовлетворяющей условию (2), можно записать системы дифференциальных и алгебраических уравнений, связанные по части переменных.
Единого названия таких систем в настоящее время не существует. В литературе используют термины: алгебро-дифференциальные системы (АДС), дескрипторные системы, системы леонтьевского типа (см., например, соответственно [1 - 3]). В зарубежной литературе утвердилось название: дифференциально-алгебраические уравнения (ДАУ)(см., например, [4]). В теории бесконечномерных задач уравнения с необратимым оператором при производной искомой функции называют уравнениями соболевского типа (см., например, [5]). Интерес к ДАУ стимулируется проблемами математического моделирования в прикладных областях, в частности, в теориях электронных схем и электрических цепей, механике и теплотехнике (см., например, [5 - 8] и приводимую там библиографию). В данной работе проведено обоснование декомпозиции одной разностной схемы для ДАУ, и эта расчетная схема применена к математической модели главного тракта прямоточного парового котла.
1. Теоремы существования
Приведем ряд необходимых для дальнейших рассуждений сведений. В работе используются нормы n-мерного вектора v = (vi V2 ••• vn )т, и х п)-матрицы V = (Vj, i =
1, 2, ■ ■ ■ , n, j = 1, 2, ■ ■ ■ , ^), вычисляемые по правилам
n
|| v У = max{|vj|, i = 1, 2, ••• , n}, || V || = max^ |vj |, i = 1, 2, ■ ■ ■ , ^}.
j=i
Под символами ||v(w)||, ||V(w)|| понимаются нормы вектор-функции v(w) и матрицы V(w), вычисленные в точке w € D С Rq. Включения v(w), V(w) € Cl(D) означают, что все частные производные элементов вектор-функции v(w) или матрицы V(w) имеют непрерывные частные производные порядка до l включительно по всем компонентам вектора w в любой точке области D. Непрерывности соответствуют включения: v(w), V(w) € C(D). Если v(w), V(w) € C(D), то их нормы также непрерывны в D.
ДАУ обладают сложной внутренней структурой. Мерой сложности является целочисленная величина, называемая индексом.
Определение 1. Пусть для системы
Л1(ж) := A(x,t, v*)x + B(x,t, v*)=0, t € T, (4)
где v* € N, существует дифференциальный оператор
l
ЛДу) :=£Lj(t,y, ■■■ ,y(l+1))(d/dt)j, j=0
со где Lj(t, v, ■ ■ ■ ,vi+i) — (n х п)-матрицы из C(T х Rn(l+2)), со свойством
My) о Л1(у) = A(y, t)y + B(y, t) Vy = y(t) € Cl+1(T),
причем det A(a(v*),a) = 0. Минимально возможное l называется индексом, системы (4) в окрестности точки (a(v*), a, v*).
Вычисление индекса является непростой задачей. Для этого нам потребуются следующие определения и утверждения.
Определение 2. (см., например [8]). Матрица, обозначаемая в дальнейшем как S-, на-
зывается полуобратной к матрице S, если она удовлетворяет уравнению: SS-S = S.
Лемма 1. [8]. Полуобратная матрица определена для произвольной матрицы S. Если выполнено условие Кронекера-Капелли: rank S = rank(S|u), то система уравнений Sy = и разрешима, и все ее решения описываются формулой: y = S-u + [E — S-S]C, где E—единичная матрица подходящей размерности, C—произвольный вектор.
Лемма 2. Пусть пучок постоянных матриц A(A) = AA + B регулярен: det A(A) ф 0. Тогда rankA > degdet A(A), где deg—символ степени многочлена.
Доказательство. Если пучок матриц A(A) регулярен, то существуют постоянные неособенные матрицы P, Q , со свойством
PA<A)Q = A(E“ ж) + (о
E
n-d
(5)
0
где Nк = 0, J—некоторый блок размерности (d х d) [9, c.354]. Отсюда
degdet A(A) = deg[det(AEd + J) det(AN + En-d)] = deg det(AEd + J) ■ 1 = d.
Следовательно, d < rankA, d = rankA тогда и только тогда, когда N = 0. □
Определение 3. Говорят, что пучок квадратных матриц AA(w) + B(w), w € D С Rq, где A— скалярный параметр (в общем случае комплексный), удовлетворяет критерию «ранг-степенью в области D, если выполнены условия
1. max rankA(w) = r, w € D; 2. det[AA(w) + B(w)] = Ara0(w) +-----, a0(w) = 0 Vw € D.
Следствие 1. Для пучка матриц AA(w) + B(w), w € D, удовлетворяющего критерию «ранг-степенью справедливо равенство: rank A(w) = const = r Vw € D.
Лемма 3. [1] Матричный пучок P(A; t) = A ^A10w)^ + , w € D, где блок A1(w)
имеет полный ранг для любого w € D, удовлетворяет критерию «ранг-степенью тогда и только тогда, когда
detP(A; w) = Ara0(w) + ..., a0(t) = det ^= 0 Vw € D. (6)
Теорема 1. Пусть для задачи (1), (2) выполнены условия:
1. A(v,t, v*), B(v,t, v*) € Cm(T х Rn), v* € N, m > 2;
2. rankA(v, t, v*) = r = const в окрестности точки ( = (a(v*), a);
3. rankA(Z, v*) = rank A(Z, v*)| b , b = —B(Z,v*);
4. rankA(Z, v*) = degdet[AA(Z, v*) + 0(c)], где
d
0(c) = B(Z,v*), B(v,t, v*) = — [B(v,t, v*) + A(v,t,v*)c] , c : A((,v*)c = —b.
Тогда:
1. индекс системы равен 1;
2. существует отрезок Т = [а, а + в] С Т, на котором определено единственное решение ж(£, V*) € Ст(Т£) задачи (1),(2).
Следствие 2. Пусть
1. A(v,t,v), B(v, t, v) € Cm(V), m > 2;
2. rankA(v, t, v) = r = const в окрестности точки (a(v*), a, v*);
3. для любого v* € N выполнены условия 3,4 теоремы 1.
Тогда существует отрезок T С T, на котором определено решение x*(t, v) € Cm(T£ х N)
задачи (1),(2).
Следствие 3. Пусть
1. A(v,t, v), B(v,t, v) € Cm(V), m > 2;
2. при v = v* на отрезке T определено решение задачи (1), (3) x(t, v*);
3. rankA(z(t)) = degdet[AA(z(t)) + 0(z(t))] = r = const, t Є T, где z(t) = (x(t, v*), t, v*),
Тогда существует окрестность
U = {w : ||w — a(v*)|| < g, rankA(w, a, v*) = rank [A(w, a, v*)|B(w, a, v*)]}
такая, что для любого b Є U на T определено решение системы (1) с начальным данным x(0, v*) = b.
Теорема 1 и следствия к ней являются частным случаем утверждений из [10]. Замечание 1. В теореме 1, согласно лемме 1, c = S-b + [E — S-S]C, S = A(a(a), a, v*). Замечание 2. Если старший коэффициент многочлена в следствии 3
обращается в нуль в точке : аг(¿*) = 0, то эта точка является особой. В ней могут
ответвляться другие решения системы (1) или другие решения могут иметь разрывы.
Пример 1. Рассмотрим две системы
где Ь Є Т = [— 1, 1]. Обе системы имеют тривиальное решение. Формально общие решения можно записать так
ж(Ь, с) = сі (в* в*)Т /(б* - 1), у(Ь, с) = ^і(Ь)сі + ^2(Ь)с2,
где сі, С2 Є И,, ^і(Ь) = {0, Ь Є Ті; #(Ь), Ь Є Т2}, ^(Ь) = {#(£), Ь Є Ті; 0, Ь Є Т2}, $(Ь) = (і2 Ь3) , Ті = [—1, 0], Т2 = (0,1], Т—символ транспонирования.
У системы 1) все другие решения имеют разрыв в точке Ь = 0. У системы 2) в этой точке ответвляются ненулевые решения, причем ранг матрицы А(Ь) постоянен и не несет информации о наличии на Т особых точек. Здесь £і(Ь; А) = (—б* + 1)А + 1, ^(Ь; А) = — ¿А + 1, аі(0) =0 в обоих многочленах. Эти нули совпадают с особыми точкам ДАУ.
К сожалению, ДАУ, у которых пучок матриц Якоби на решении удовлетворяет критерию «:ранг-степень>, не исчерпывают всех систем индекса 1.
Пример 2. Рассмотрим тестовую задачу вида (1), (3)
окрестности начальной точки. Тем не менее, простое вычисление показывает, что в определении 1 можно принять
д
0(z(t)) = B(z(t)), B(v,t, v*) = — [B(v,t,v*) + A(v,t, v*)c], c = x(t, v*).
{(t; A) = det[AA(z(t)) + 0(z(t))] = ar(t)Ar + ■ ■ ■ + ao(t)
(7)
где (^1(0) х2(0))Т = (0 0)Т . Матрица А(ж,£, V) здесь меняет ранг в сколь угодно малой
У1 ) Vy Є C2(T).
В различных областях приложений применяют следующий прием: ДАУ регуляризируют, дифференцируя алгебраические связи. В нашем случае это приводит в тупик. После дифференцирования в системе (7) второго уравнения X2 =0 мы получаем задачу, у которой матрица перед производной вырождается на решении, и непонятно, что делать дальше.
Достаточно часто в приложениях дифференциальные и алгебраические уравнения разделены, и система (1) имеет вид
A(t,x,v)x + B(t,x,v) = (Al(t0X,V^ x + (Bllt’X’V)) =0, t Є T = [а,в]‘ (8)
Для системы (8) справедливо утверждение о разрешимости ДАУ на всем отрезке T. Теорема 2. Пусть для задачи (8), (2) выполнены условия:
1. A(v,t, v), B(v,t, v) Є Cm(V), m > 2;
2. max{rankA(v, t, v), (v,t,v) Є V} = r;
3. rankA(v, t, v) = rank A(v, t, v)| b(v) , b(v) = — B(v,t,v);
4. старший коэффициент многочлена
dB(v, t, v)
det
AA(v, t, v) +
= ar (v, t, v)Ar +
dv
удовлетворяет условию |ar(v,t, v)| > co > 0, V(v,t, v) Є V;
5. ||Ai(v,t, v)|| + ||dB2(v,t, v)/dv^ < Ki, ||B(v,t,v)|| < k2 + кз|М|,
||dB2(v, t, v)/dt|| < k4 + k5IIv|| V(v, t, v) Є V, K = const, i = 1, 5.
Тогда существует единственное решение x(t, v) Є Cm(T x N) задачи (8),(2). Доказательство проводится тем же способом, как в теореме 1 из [10].
2. Численные методы
Обоснование численных методов для решения задачи (1), (2) связано с со значительными трудностями. На отрезке T = [а, в] введем сетку ti = а + ih, i = 0,1, ■ ■ ■ , M — 1,
h = (в — a)/M, M - число узлов сетки. Запишем неявную схему Эйлера в двух вариантах
A(Y+i) Xi+1h— Xi + B (Y+i) = 0, (9)
А(Гі) Xi+1h— Xi + B(Yi+i) = 0, (10)
где Y = (ti, Xi, v*), Xo = a(v*), и при каждом i нужно один раз решить нелинейную систему.
Пример 3. Рассмотрим тестовую задачу вида (1), (3)
X1 Х2 Vх Л + (X1 — Х2 — 2еЛ =0 t Є T =[01] (11)
(xi — Ж2)Ж1 (xi — X2)X2j \X2) V (xi — X2)2 ) , [ ( )
где (жі(0) Ж2(0))Т = (1 1)Т . Применим теорему 1. Имеем
А« ><= = (І І) с = -в(< )=(0) •« =(2 V« € к,
гапк А(() = degdet[ЛA(Z) + 0(с)] = degdet
. , 1 1\ /1 + к 1 — к
Л[ 0 0І+1 2 —2
= deg(—4Л — 4) = 1.
Согласно теореме 1, начальная задача для системы (9) имеет решение. Методом исключения неизвестных можно найти решение (жі(і) Ж2(і))Т = (е* е*)Т . Применим схемы (9), (10) для решения системы (9). Нелинейные системы будем решать методом Ньютона. За начальные данные при вычислении Жі+і будем брать вектор Жі [11]. Легко видеть, что для системы (11) в методе (9) матрица Якоби вырожденная для любой итерации, и метод Ньютона не применим.
Проведем анализ схемы (10). Она гораздо сложнее чем схема (9). Упростим ее. Разлагая в ряд Тэйлора вектор-функцию В(Уі+і) и матрицу А(Уі+і) в точке Жі приведем (9) к виду
[Аі + Н(Ві + Аі](жі+і — Жі) = —^вії (12)
где А = ■(!•), В(І-),
~ дВ(-и,^ V) ~ д [А(-и,і, V )Рі]
Ві = д^ |(^,і,^)=Уі, Аі = — |(^,і,^)=Уі, Рі = (жі — Жі-і)/Н, і ^ 1-
Если і = 0, то принимаем Ро = с, где с—вектор из теоремы 1.
Теорема 3.
1. Пусть выполнены условия теоремы 1;
2. на отрезке Тє = [а, а + е] задана сетка ¿і = а + іЛ, і = 0,1, ■ ■ ■ , М — 1, Н = е/М.
Тогда, начиная с некоторого М > Мо, система (12) разрешима для всех і относительно Жі+і, и справедлива оценка ||жі — ж(^)|| = О(Н) равномерно по і.
Метод ниже не используется, поэтому мы опускаем доказательство.
Лемма 4. (см., например [12]) Пусть для неотрицательных чисел щ, і = 2, 3,... выполнено неравенство
і-і
щі < с + т ^ -щ, щі < с, і=і
где с, т - неотрицательные постоянные. Тогда справедлива оценка (разностный аналог
леммы Гронуолла-Беллмана): щ < с(1 + m)i-1, і = 1, 2,....
Для схемы (10) эти трудности отсутствуют. Проведем следующие преобразования. Продифференцируем второе уравнение в системе (8) и выпишем разностную аппроксимацию этого выражения
5 Жі+і — Жі _ п 5 ^
В2,і н + ^і = 0, В2,і = дЖ |(^,*,^)=уі, ^і = д£ |(^,*,^)=уі•
Пусть выполнены условия теоремы 2. Из леммы 3 вытекает, что det £і = det (■■і’М = 0 Уі,
\В2,і/
и матрица В?2,і имеет полный ранг для любого і. Тогда можно произвести разбиение В?2,і =
B22,i) где det B22,i = 0 для любого i. Согласно разбиению матрицы B2,i, произведем разбиение: x = (yT zT)T. Тогда из схемы (10) следует
у! yi + 1 yi I /" \ П Zi + 1 о — 1 f D yi + 1 yi I 1 /-I
Ai h + B1(yi,zi,v*) = 0J h = —B22,i[B21,i h + Gi]’ ( )
Здесь аргумент при вычислении функции B сдвинут на один шаг. Матрица Ai = Ац,і — A12,iB—2^21^ в (13) получена применением формулы Шура [9] к матрице Si, следовательно, det Ai = 0 Vi. Итак, мы произвели декомпозицию схемы (10), выписав схему (13). Здесь нам нужно обратить две матрицы Ai, B22,i. Мы получили существенное уменьшение объема вычислений. Например, если размерности матриц Ai, B22,i равны, то количество операций уменьшается в 4 раза по сравнению с применением исходной схемы (10). К тому же нам не нужно решать нелинейные уравнения. Это очень важно при проведении вычислений в режиме реального времени.
Если мы подставим решение задачи (8), (2) в первое уравнение из (13), вычтем из полученного выражения уравнение (13) и перейдем к нормам, то в условиях теоремы 2 получим
i— 1
выражение Wi < а + wj, w = ||z(ti) — Zi||, а = O(h), к = const. Из леммы 4 следует,
j=1
что справедлива оценка Wi < ек(в—1а)а.
Замечание 3. Если система (8) жесткая, то в (13) можно сдвиг аргумента производить следующим образом: слагаемое B1 (yi,Zi,v*) заменить на слагаемое В.(уі+1, Zi, v*), и мы получим неявную схему. На сходимость это не влияет.
3. Модель основного тракта прямоточного парового котла
В технологической установке, называемой паровым котлом, вода в сети трубопроводов (схематично изображенных на рис. 1) последовательно нагревается в устройствах, называемых теплообменниками (ТО), по ходу от насоса к турбине газами и лучевым теплом, получаемыми от сгорания топлива в топке. Газы идут в в противоположном направлении по отношению к ходу воды.
Я
*-------•
1 р1 2 р2 3
Рис. 1. Принципиальная схема главного тракта прямоточного котла
В ТО 1 происходит нагревание воды до температуры кипения, в ТО 2 (зона кипения) вода превращается в пар, а в ТО 3 (участок перегрева) пар доводится до нужной температуры и давления. По ветке 4 в ТО 3 подается вода из ТО 1 для регулирования температуры пара на выходе из котла. Содержание тепла в среде характеризуется величиной, называемой энтальпией и обозначаемой ниже буквой I. В простейшем случае I = ct, где с— теплоемкость среды, а t— ее температура. Система уравнений, описывающая ТО имеет вид
дГ- —
V' (1j — ) = — t(pi> )]; (14)
. ___________________________________________________
Mm,jcm"dTT = (tg,j — ) — [^j — t(pi, )], (15)
pg cg jgj = Dg cg (tS,j ¿gbxj ) ag,j Hg,j (tg,j )• (16)
dtg dr
где j = 1,2,3— номер ТО, i = 1,2, , tgbxj — входные в ТО энтальпия и температу-
ра газа, Dbj, Dg — расходы теплоносителя (пара, воды или пароводяной смеси) и газа через ТО, Ij — энтальпия на выходе из ТО, Ij — средняя энтальпия, Ij = 7/j + (1 — Y)Ibx, 7 € [0,1]—коэффициент усреднения t(p, I)— функция, связывающая температуру теплоносителя с энтальпией и давлением в ТО; р = р(р, I)— функция плотности теплоносителя в ТО; - температура стенки ТО, tgj —температура газа в ТО, Dg - расход газа, cm, cg — теплоемкости газа и металла, Vj, Vgj, Hbj, — объемы ТО и площади тепловоспринима-
ющих поверхностей по воде и газу, abj, agj — коэффициенты теплоотдачи по теплоносителю и газу, qLj — лучевое тепло, воспринимаемое ТО.
В основу моделирования топки положен закон Стефана - Больцмана.
dT
Vgpg= -Dgcg(Tg - rCK) - + Qug1, (17)
Mzc*-^ =---------(Tz - Ti) + Ql, Ql = aiH([z(Tg)]4 - T4), (18)
dTz H
^ =-------(
dr pz
где Tg = ig, + rCK, Tz = + rCK, T1 = t(p2,I2) + rCK, rCK = 273.15 - соответственно тем-
пературы газа на выходе из топки, тепловоспринимающего слоя и теплоносителя в градусах Кельвина, z = cTg - температура факела, c = const = 1,1, Ql — лучевое тепло, выделяемое топкой, Qug —тепло выделяемое от сгорания топлива (угля), Ql = + ^¿,2 + qL,3- Система уравнений, описывающих математическую модель ТО или топки, выражает законы
сохранения. Количество тепла, поступающее из газом и лучевым теплом в теплообменник Dgcg (tgj — ) + qL,j , в стационарном состоянии равно количеству тепла, уносимого теп-
лоносителем Dbj (Ij — /5x j)• Для топки количество тепла от сгорания топлива равно теплу, уносимому с газом и лучевым теплом: Qug = Dgcg(Tg — rCK) + Ql.
Сеть трубопроводов моделируется гидравлической цепью (ГЦ). Под гидравлической цепью понимается плоский граф (для примера см. рис.1), для которого выполнены первый и второй законы Кирхгофа: а) количество среды втекающей в вершину графа, равно количеству вытекающей; б) сумма перепадов давлений в замкнутом контуре равна нулю. ГЦ из m узлов и r ветвей ставится в соответствие (m х r) —матрица соединения ветвей и границ A, состоящая из нулей и единиц [13]. Уравнения, вытекающие из этих законов, дополняются соотношениями, связывающими давления в вершинах графа с расходами по ветвям. Используя результаты из [7], нестационарную модель ГЦ запишем в виде ДАУ
(Я 0\ /D(t)WSo AT\ (D(t)\ + (S|D(i)|D(t)\ = (H(i) + ATP,(i)\ . _ [0 ^ (|9)
1.0 0J UW4 A 0) U‘)J4 0 M Q(t) )'1 € [0 <19>
где (AT Aj") = A, D(t) = (di(t) d2(t) , ••• , dr(t))T - вектор-функция расходов среды по ветвям, P(t) = (pi(t) P2(t), ••• ,P^(î))T - вектор-функция давлений в узлах,
P*(t) = (Pi,*(t) P2,*(t), ••• ,Pv,*(t))T - вектор-функция известных давлений, ^ + V = m, Я = diag{^i,^2, ••• , £г} > 0 - параметры инерции ГЦ, S0 = diag{si;0, s2,0 ••• ,sr,0} > 0, S = diag{si,s2, ••• , sr} > 0 - параметры сопротивлений ветвей ГЦ, |D(t)|D(t) = (di(t)|di(t)|, d2(t)|d2 (t)|, ••• ,dq (t)|dq (t)|)T, H (t) = (hi(t) h2(t), ••• , h (t))T, Q(t) =
(qi(t) dq(t), ■ ■ ■ , qr(t))T - вектора напоров и притоков в ГЦ. Индекс системы (19) глобально равен 2. В определения 1 можно принять
Ло
Ег 0\ d ( 0 0
0 0у dt \—A
R-i 0\ + _d ^0 0
0 0^ + dt V0 Е
/ Ег 0
l^(D) —AR-iAT
о
где ёе1 АЛ-1АТ > 0, Ф(^) = — [АЛ-1£0 + д(А£|^|^)/д^]. При нулевых Н(¿), Р*(£) и ф(£) система (19) имеет нулевое решение, которое является контрактивным. Рассмотрим скалярное произведение вектор-функции (^Т(¿), —РТ(¿))Т с системой. Получим
й
-(£(*),£ т(*)) = -5о(Я(*),£>т(*)) - 5(|Я(*)|Я(*)|,£>т(*)) < 0.
Следовательно, ||^(£)|| ^ 0, í ^ то. В силу алгебраических связей ||Р(і)|| ^ 0, і ^ то.
Замечание 4. Если в (19) положить Я = 0, то мы получим стационарную ГЦ, которая исчерпывающе изучена в [13].
Используя формулы (14) - (18), выпишем ДАУ, описывающее модель котельной установки, принципиальная схема которой изображена на рис. 1. Параметры соответствуют реальной установке (прямоточному котлу ПК-24). Элементами управления здесь служат параметры и 1ьх — теплосодержание воды на входе в котел, сопротивления участков
$3, $4, регулирующие количество воды, поступающее в котел и температуру пара на входе в турбину.
Уравнения ГЦ: <
Уравнения TO 1: Уравнения TO 2:
Уравнения TO 3:
0 = жз — ж4 — Жб, q(X2, Ж10) = Ж4 — Ж5 + Жб,
£зЖ 3 = Ж1 + 80,зЖз + S3 + К1 /0*(545 — X13(s)ds |Жз|Жз — h.1 (t),
^4^Ж 4 = —Ж1 + Ж2 + S0,4X4 + 54|Ж4|Ж4,
£5Ж 5 = —Ж2 + ¿0,5Ж5 + S5 |Ж5 |Ж5 — h2(t)
^бЖ б = —Ж1 + Ж2 + 50,бЖ4 + 8б|Жб|Жб;
—а11Ж 7 = а12Ж3 (ж7 — 4x0 — а13(ж8 — а14 Ж7), а15^Ж 8 = «1б(Ж9 — Ж8) + Й1з(Ж8 — «14 Ж7),
&17Ж 9 = «18 (Ж9 — Ж12) + «1б (Ж9 — Ж8)
—Й21^Ж 10 = Й22Ж4(Ж10 — Ж7) — Й2з(Жц — Й24Ж10),
—а25^Ж 11 = —Й2б(Ж12 — Жц) — 023(Ж11 — Й24Ж10) — ^ф(Ж1б, Ж17), 027Ж 12 = «28 (Ж12 — Ж15) + 02б(Ж12 — Жц) а31Ж 1з = а32ж5(ж13 — Tsm) — а33(ж14 — a34x13^
—аз5Ж 14 = —азб(ж15 — Ж14) + азз (Ж14 — 034X13) — ^2^(ж1б, Ж17), —а37Ж 15 = а38 (ж15 — ^з(ж1б — rCK)) + а3б(ж15 — x14^
Tsm — (Ж6Ж7 + ж10ж4)/ж5;
16 — Й42(Ж16 - rCK) + Q(X16, Ж17) - Qugl,
Уравнения топки: ^ а45ж 14 — а46(ж17 — ж11) — Q(x16, ж17),
Q(X16,ж17) — Й43 [(а44ж16)4 — ж^] .
Предполагается, что вода в ТО 2 находится на линии насыщения. Приток д(ж2,жю) —
— V2dp/dt появляется в результате учета изменения плотности воды в ТО 2 при возмущениях. В стационарном состоянии д(ж2,жю) — 0. Если предполагать что гидравлические процессы протекают существенно быстрее, чем тепловые (в (19) R — 0), то система вблизи стационарного существования имеет индекс 1. В обратном случае индекс системы равен 2. Укажем соответствие между физическими параметрами и переменными системы. Здесь
(ж1, ж2, ж3, ж4, ж5, ж6, ж7, же, ж10, ж11ж12, ж13, ж14, ж15, ж16, ж17) —
(Р1,Р2, db ^ ^ d4, #ъ ¿й,Ъ ^ ^ ^ #3, Tg, Tz)-
Таблица
Основные параметры, их обозначение и значения
Обзначение Численная величина Обзначение Численная величина
«11 500 «21 1000
«12 0,5 «22 0,5
«13 102.24 «23 100
«14 1/4,5 «24 0,1
«15 50 000 «25 100 000
«16 690 «26 700
«17 100 «27 100
«18 132,78 «28 140
«31 1000 «41 3700
«32 0,5 «42 139
«33 100 «43 0, 265110222 х 10-7
«34 1/7,5 Й44 1,1
«35 100 000 «45 825
«36 700 «46 1125/7,8
«37 100 ^(¿) 200
«38 140 Н2 (¿) 140
«3 0,002 «0,3 0
«4 0,002 «0,4 0
«5 0,002 «0,5 0
«6 0,1 «0,6 0
Х1(0) 167 жю(0) 2600
Х2(0) 154 жи(0) 590
Х3(0) 100 Ж12 (0) 600
Ж4(0) 95 Ж13(0) 4000
Жб(0) 5 Ж15 (0) 650
Ж7(0) 1193 Ж16(0) 1379,266638
Ж8(0) 400 Ж17(0) 1202,612331
Жд(0) 420
Проведены расчеты при К = 0 с использованием декомпозиции и без на временном отрезке [0,1200] (время измеряется в секундах) с шагом интегрирования Н = 0, 05, коэффициент усреднения 7 = 0, 5. Результат расчетов с декомпозицией изображен на графике ниже (рис. 2), который иллюстрирует поведение температуры пара на выходе из котла (переменная Ж1э/Ср, ср - теплоемкость пара). Начальные данные намеренно взяты далеко от стационарных. Температура на выходе из котла должна быть равна 545 °С. В масштабах графика различия (с декомпозицией и без) не видны. Точное сравнение показало, что максимальная величина расхождения расчетов с декомпозицией и без нее равна 0,2762.
Расчеты при К > 0 проведены по следующей схеме: индекс системы понижался до 1, а затем применялась схема с декомпозицией. Понижение индекса проводилось по методу из [14]. Результаты вычислений для нестационарного случая ГЦ отличаются незначительно от случая К = 0. При всех вариантах принято = 294074, 55, /^ = 1000.
Работа проводилась при финансовой поддержке РФФИ, гранты №09-08-00201-а, 11-01-93005-Вьет-а
640 620 600 580 560 540 520 500 480 460 440 420 400
0 100 200 300 400 500 600 700 800 900 1 ООО 1 100 1 200
Рис. 2. Поведение температуры пара на выходе из котла Литература
1. Бояринцев, Ю.Е. Алгебро-дифференциальные системы. Методы численного решения и исследования /Ю.Е. Бояринцев , В.Ф. Чистяков. - Новосибирск: Наука, 1998.
2. Зубова, С.П. Решение задачи управления для линейной моделей дескрипторной системы с прямоугольно-матричными коэффициентами / С.П. Зубова // Матем. заметки. - 2010.
- T. 88, № 6. - С. 885 - 896.
3. Келлер, А.В. Численное исследование задач оптимального управления для моделей леонтьевского типа: дис. ... д-ра физ.-мат. наук /А.В. Келлер. - Челябинск, 2011.
4. Хайрер, Э. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи / Э. Хайрер, Г. Ваннер. - М.: Мир, 1999.
5. Свиридюк, Г.А. Задача Коши для линейного сингулярного операторного уравнения типа Соболева / Г.А. Свиридюк //Дифференц. уравнения. - 1987. - T.23, № 126. -
С. 2169 - 2171.
6. Ушаков, Е.И. Статическая устойчивость электрических систем / Е.И. Ушаков. - Новосибирск: Наука, 1988.
7. Балышев, О.А. Анализ переходных и стационарных процессов в трубопроводных системах (теоретические и экспериментальные аспекты) /О.А. Балышев, Э.А. Таиров. -Новосибирск: Наука, 1998.
8. Бояринцев, Ю.Е. Регулярные и сингулярные системы линейных обыкновенных дифференциальных уравнений /Ю.Е. Бояринцев . - Новосибирск: Наука, 1980.
9. Гантмахер, Ф.Р. Теория матриц / Ф.Р. Гантмахер. - M: Наука, 1967.
10. Чистякова, Е.В. О разрешимости вырожденных систем квазилинейных интегро-дифференциальных уравнений общего вида /Е.В. Чистякова, В.Ф. Чистяков // Вычислительные технологии. - 2001. - Т. 16, № 5.-C. 100-114.
11. Калиткин, Н.Н. Численные методы /Н.Н. Калиткин . - М.: Наука, 1978.
12. Апарцин, А.С. Неклассические уравнения Вольтерра I рода: теория и численные методы /А.С. Апарцин. - Новосибирск: Наука, 1999.
13. Меренков, А.П. Теория гидравлических цепей парогенераторов / А.П. Меренков, В.Я. Хасилев. - М.: Наука, 1985.
14. Булатов, М.В. О преобразовании алгебро-дифференциальных систем уравнений / М.В. Булатов // Журн. вычисл. математики и мат. физики. - 1994. - Т. 34, № 3. - C. 360 - 372.
Виктор Филимонович Чистяков, доктор физико-математических наук, главный научный сотрудник, Институт динамики систем и теории управления СО РАН (Иркутск, Российская Федерация), chist@icc.ru.
Эмир Асгадович Таиров, доктор технических наук, заведующий лабораторией, Институт систем энергетики им. Л.А. Мелентьева СО РАН (Иркутск, Российская Федерация), tairov@isem.sei.irk.ru.
Елена Викторовна Чистякова, кандидат физико-математических наук, научный сотрудник, Институт динамики систем и теории управления СО РАН (Иркутск, Российская Федерация), elena.chistyakova@icc.ru.
Анатолий Алексеевич Левин, кандидат технических наук, старший научный сотрудник, Институт систем энергетики им. Л.А. Мелентьева СО РАН (Иркутск, Российская Федерация), levin@isem.sei.irk.ru.
On Decomposition of Difference Schemes for Numerical Solution of Differential Algebraic Equations
V.F. Chistyakov, Institute for System Dynamics and Control Theory of Siberian Branch of Russian Academy of Sciences (Irkutsk, Russian Federation),
E.A. Tairov, Melentiev Energy Systems Institute of Siberian Branch of the Russian Academy of Sciences (Irkutsk, Russian Federation),
E. V. Chistyakova, Institute for System Dynamics and Control Theory of Siberian Branch of Russian Academy of Sciences (Irkutsk, Russian Federation),
A.A. Levin, Melentiev Energy Systems Institute of Siberian Branch of the Russian Academy of Sciences (Irkutsk, Russian Federation)
We consider quasi-linear systems of ordinary differential equations (ODE) with identically singular matrix multiplying the derivative of the desired vector-function and difference scheme for their numerical solution. We discuss conditions that make it possible to solve algebraic (finite-dimensional) equations at each step of numerical process and substitute the solutions obtained into the dynamics equations. Results of numerical solution of ODE systems modeling direct-flow boiler unit are given.
Keywords: differential-algebraic equations, index, difference schemes, mathematical models, direct-flow boiler.
References
1. Boyarintsev Yu.E., Chistyakov V.F. Algebro-differentsial’nye sistemy. Metody chislennogo resheniya i issledovaniya [Differential Algebraic Equations. Methods of Numerical Solution and Research]. Novosibirsk, Nauka, 1998.
2. Zubova S.P. Solution of the Control Problem for a Linear Descriptor System with Rectangular Coefficient Matrix. Mathematical Notes, 2010, vol. 88, no. 5-6, pp. 844 - 854.
3. Keller A.V. Chislennoe issledovanie zadach optimal’nogo upravleniya dlya modeley leont’evskogo tipa: dis...doktor fiz.-mat. nauk. Chelyabinsk, 2011.
4. Hairer E., Wanner G. Reshenie obyknovennykh differentsial’nykh uravneniy. Zhestkie i differentsial’no-algebraicheskie zadachi [Solution of Ordinary Differential Equations. Rigid and Differential-algebraic Problems]. Moscow, Mir, 1999.
5. Sviridyuk G.A. The Cauchy Problem for a Linear Singular Operator Sobolev-type Equations [Zadacha Koshi dlya lineynogo singulyarnogo operatornogo uravneniya tipa Soboleva]. Differential Equations, 1987, vol. 23, no. 126, pp. 2169 - 2171.
6. Ushakov E.I. Staticheskaya ustoychivost’ elektricheskikh sistem [Static Stability of Electrical Systems]. Novosibirsk, Nauka, 1988.
7. Balyshev O.A., Tairov E.A. Analiz perekhodnykh i statsionarnykh protsessov v truboprovodnykh sistemakh (teoreticheskie i eksperimental’nye aspekty) [Analysis of Transient and Steady-state Processes in Pipeline Systems (Theoretical and Experimental Aspects)]. Novosibirsk, Nauka, 1998.
8. Boyarintsev Yu.E. Regulyarnye i singulyarnye sistemy lineynykh obyknovennykh differentsial’nykh uravneniy [Regular and Singular Linear Systems Ordinary Differential Equations]. Novosibirsk, Nauka, 1980.
9. Gantmakher F.R. Teoriya matrits [Matrix Theory]. Moscow, Nauka, 1967.
10. Chistyakova E.V., Chistyakov V.F. On the Solvability of Degenerate Systems of Quasilinear Integro-differential Equations of the Form [O razreshimosti vyrozhdennykh sistem kvazilineynykh integro-differentsial’nykh uravneniy obshchego vida]. Computational Technologies, 2001, vol. 16, no. 5, pp. 100 - 114.
11. Kalitkin N.N. Chislennye metody [Numerical Methods]. Moscow, Nauka, 1978.
12. Apartsin, A.S. Neklassicheskie uravneniya Vol’terra I roda: teoriya i chislennye metody [Non-classical Equations of Volterra Type I: Theory and Numerical Methods]. Novosibirsk, Nauka, 1999.
13. Merenkov A.P., Khasilev V.Ya. Teoriya gidravlicheskikh tsepey parogeneratorov [Theory of Hydraulic Circuits Steam]. Moscow, Nauka, 1985.
14. Bulatov M.V. The Transformation of Differential-algebraic Systems of Equations [O preobrazovanii algebro-differentsial’nykh sistem uravneniy]. Computational Mathematics and Mathematical Physics, 1994, vol.34, no. 3, pp. 360 - 372.
Поступила в редакцию 19 декабря 2011 г.
