CC BY
52
УДК 622.691 - 50 ББК 22.161:33.362
МНОГОСЕТОЧНЫЕ БАЛАНСОВЫЕ МОДЕЛИ НЕСТАЦИОНАРНЫХ ПОТОКОВ В СЛОЖНЫХ ГАЗОТРАНСПОРТНЫХ СИСТЕМАХ
12 3
Ахметзянов А. В. , Сальников А. М. , Спиридонов С. В. ,
(Учреждение Российской академии наук Институт проблем управления РАН, Москва)
Предлагаются многосеточные варианты иерархических методов моделирования распределения потоков для оптимизации нестационарных режимов транспорта газа на различных уровнях диспетчерского управления Единой системой газоснабжения (ЕСГ) России. Для увеличения точности вычислений многоуровневых методов на каждом уровне иерархии вычислений используется многосеточная схема, основанная на релаксационной схеме Р. П. Федоренко [4, 5]. Универсальность технологии определяется легкостью построения многосеточной структуры и операторов перехода, обеспечивающих высокую точность дискретизации на грубых сетках без предварительного сглаживания и интерполяции, присущих классическим многосеточным методам [2, 3]. При конечно-объемной дискретизации универсальная многосеточная технология пригодна для объектов, описываемых уравнениями с разрывными и нелинейными коэффициентами и правыми частями.
Ключевые слова: многосеточный метод, балансовая модель, нестационарный поток, транспорт газа.
1 Атлас Валиевич Ахметзянов, кандидат технических наук (Москва, ул. Профсоюзная, д. 65, тел. (495) 334-90-30, awa@ipu.ru).
2 Антон Михайлович Сальников (salnikov@ipu.ru).
3 Сергей Владимирович Спиридонов (spirid@ipu.ru).
1. Введение
Газотранспортные системы (ГТС), образующие Единую систему газоснабжения (ЕСГ) страны, представляют собой сети газопроводов произвольной конфигурации (с закольцованными подсистемами) и предназначены для магистрального транспорта газа от месторождений к промышленным и административным центрам, включая экспорт газа в Европу. Отдельные ветви газопроводов - это несколько параллельных ветвей трубопроводов, соединенные между собой перемычками. Вдоль ветвей магистральных газопроводов могут быть расположены промежуточные компрессорные станции (КС), а также промежуточные отборы (стоки) для крупных потребителей (ГРЭС, ТЭС и др.) и притоки от источников, т. е. месторождений и подземных хранилищ газа (ПХГ).
Сложность конфигурации и большая протяженность рассматриваемого объекта математического моделирования с использованием обычных методов вычислительной математики (например, конечно-разностных) приводят к необходимости решения системы уравнений очень большой размерности. Проблема размерности может быть преодолена с использованием многоуровневой (иерархической) декомпозиции системы и соответствующего распараллеливания вычислительных процедур [6]. Наибольшая эффективность в смысле объема, времени и точности схемы распараллеливания вычислений достигается, если выбранная многоуровневая декомпозиция системы на отдельные подсистемы произведена с максимально возможным учетом особенности структуры рассматриваемых объектов моделирования, т. е. региональных ГТС и ЕСГ в целом. При таком подходе требуемая точность вычислительных алгоритмов и существенное уменьшение необходимого объема вычислений достигается при использовании универсальной многосеточной технологии для формирования системы уравнений на разных уровнях иерархии вычислений.
Предлагаемые методы ориентированы на использование современных средств параллельного программирования и кластерных вычислительных систем, в частности MPI (Message Passing Interface) и суперкомпьютеров.
2. Постановка и методы решения задачи для участка газопровода
При изотермическом (T = const) течении реального газа (p / р = zRT) нестационарное распределение давления p(x, t) вдоль участка неразветвленного газопровода в момент времени t описывается системой уравнений в частных производных
-dp/dx = d(pw)/dt + 1\w | (pw)/2d, 0 < x < l;
(1) -dp =- _L dp / d t = d(pw)/ d x, t > 0, d t c
где р = p / c2 — плотность газа, z - коэффициент сжимаемости газа, R - газовая постоянная, T - абсолютная температура, c -скорость звука в газе, X - коэффициент гидравлического сопротивления, d - диаметр гидравлического сечения трубы, l - длина участка газопровода, w - скорость газа.
В газотранспортных системах выделяют следующие три основных типа моделирования режима течения газа: X = 0 -незатухающее волновое или инерционное (преобладание сил инерции и слабое влияние трения), d(pw) / dt = 0 - безынерцио-ное (преобладание сил трения и отсутствие влияния сил инерции) и затухающее волновое (общий случай с учетом сил инерции и трения).
Для задач диспетчерского управления ГТС и ЕСГ в целом характерным является режим с преобладанием влияния сил трения. Это обусловлено геометрическими (большая протяженность участков с диаметрами труб в несколько метров) и технологическими (большое гидравлическое сопротивление и большой расход газа) параметрами. Поэтому (1) можно представить в виде
-дР/дx=Сс21 q | q/d, 0<x<I;
(2.1)
-дP/дt=2с2 дд/дx, t >0,
д2P 2С^|дРл , л (2.2) —=———,0<х<1, t >0,
дх2 d дt, ’ ’
где Р = р2 и q = pw — массовый расход через единичную поверхность сечения трубы. Это допустимо, поскольку волновой характер движения газа быстро затухает, а исходная система (1) заменяется эквивалентным нелинейным уравнением параболического типа.
Для преобразованной системы (2.1) должны быть заданы начально-краевые условия (любое, но только одно из ниже указанных)
Р(0, t) = ^лев, p(l, t) = Рпр ^), Р( х,0) = Р (x);
( ) Р(0, t) = PЛев, q(l,t) = qпр (t), P(x,0) = Рй(x),
где P0(x) - распределение квадрата давления вдоль участка в начальный момент времени t = 0, Pлев(t), Pnр(t), qлев(t), qпр(t) -заданные функции изменения давления и расхода на границах участка в рассматриваемом интервале времени [0, 7]
Для простоты и наглядности будем предполагать, что начальные условия соответствуют стационарному потоку, т. е. для t < 0 dP(x,t) / dt = 0. В этом случае уравнение (2.1) при начальнокраевых условиях, соответствующих второй строке (3), приобретает вид др2 / дх = -Сс2д2 / d, следовательно,
(4) Р(х,0) = Р0 (х) = д/Рлев - Сс2д2рх / dP, 0 < х < I.
Двухточечные начально-краевые задачи для нелинейного параболического уравнения (2.2) при любом традиционном варианте начальных и граничных условий являются корректно поставленными, т. е. непрерывным изменениям граничных условий (3) х = 0 и х = I соответствуют непрерывные изменения решения (2.2).
Явная и неявная сеточная аппроксимация уравнения (схема вычислений) (2.2) по времени имеет вид
р («+1) - р(я) д2р|-я+1)
(5) а--------:------= „ 2 ,
А/ дх
р («+1) - р(я) д2р(«)
(6) а--------------= . ,
А/ дх
где о = 2Х^ / d.
Конечно-объемная (интегрально-балансная) аппроксимации
(5) и (6) по пространственной координате х выражаются матричными уравнениями
(7) — ср(я+1) - — Ср(я) + Вр(я+1) = 0 « (— С + В)р(я+1) = — Ср(я),
А/ А/ А/ А/
(8) — ср (я+1) —- Ср(я) + Вр(я) = 0 « — Ср (я+1) = (-- С - В) р(я),
А/ А/ А/ А/
где матрицы В и С - матрицы, соответствующие аппроксимации
р 2С | q | др Е
членов -dгv^rad р и------------------, соответственно. Если положить,
d д /
что о = 2^(я)| / d, где распределение массового расхода газа q(я) = q(x, /)|/ = яА/, матричные уравнения (7) и (8) также становятся линейными системами алгебраических уравнений Ар(я+1) = р, где А = С / А/, р = (с/А/ + В)р(я) и А = С / А/ - В, р = Ср(я) / А/, соответственно.
Для решения уравнений конечно-объемных аппроксимаций (7) и (8) на я-м временном слое и определения q(я) = q(x, /)|г = яА/ , воспользуемся следующими схемами итерационных вычислений. Сначала на линейный участок газопровода О = [0, I] нанесем точки с координатами хк : 0 = х0 < ... < хк < ... < хк = I, т. е. построим ячейки Ок = (хк, хк + !), к = 1, К -1. Левую и правую половины ячейки Ок представим открытыми интервалами О к = (хк, (хк + хк+1) / 2) и О к = ((хк + хк + 1) / 2, хк), затем построим конечный объем Ок = О к _ 1 и Ок вокруг каждого внутреннего узла. Конечные объемы, соответствующие граничным узлам х! и хк определим как О1 = и ОК = О.'Кс-1, уравнения (5) представим в общем виде - divgradр + ар = ар, где р = р(я+:), р = р(я),
о = 2Я|д(и)| / (й • АО, и проинтегрируем его по каждому конечному объему 0,к, в результате получим следующую систему уравнений
+^Ц3Р0 + Р) = ^_(3Ро + Р)
А „ „ „ Л
Р1 - Р0 йР
к1 йх х=х0 =0 0
V
Р - Р
Гк+1 к
к Рк-1
Л
+
ой.
(3Рк + Рк+1) =
к-1' Йк-1°
+
ой,
к-1
(3Рк + Рк-1) +
ко
(3Рк + Рк-1) + ^ (3Рк + Рк+1),
8
йР
йх
Р - Р
ГК гК-1
к
К-1
°-1(3Рк + Рк-,) = ^(ЭРк- + Рк-1).
8
8
где кк = хк + 1 - хк, Р - ранее рассчитанное значение, к = 1, К -1.
После перегруппировки подобных членов получим следующее матричное уравнение:
АР = (В + С) Р = С с матрицей А = (ОкУ)к ]=1к = (*ку + ск;)к,^, содержащей следующие ненулевые элементы: а11 = 1/к1 + 3ок1 /8 , а12 = -1/к1 + ок1 /8 , аКК = 1,
а
к (к-1)
= -1/кк-1 + 3^кк-1/8, ак (к+1) =-1/кк + 3скк/8,
акк = 1/кк-1 +1/кк + 3°(кк-1 + й, )/8,
а
К ( К-1)
= -1/кК-1 +окК-1 /8, аКК = 1/кК-1 + 3окК-1 /8 .
А вектор правой части С = £ £2, ..., gk) содержит элементы £1 = к°(3Роо + Р)/8,
gk = кк-о(3Рк + Рк-1) / 8 + кко(3Рк + Рк+1) / 8, "к е {1,!,К -1},
£К = кК-1°К-1/2(3РК-1 + РК )/8.
Отсюда следует, что ненулевые элементы матриц В и С определяются формулами
Ьи = 1/к1,Ь12 = -1/к1 и с11 = 3ок1 /8, с12 = ок1 /8 ;
к
8
к
8
8
+
х=х^=А
К
ьк (к -і) = -1/К-і, ьк (к+і) = 3о\/8; и
Ск(к-1) = 3оКк-і /8, Ск(к+1) = -1/Кк ;
ькк =1/Кк-і +1/К и скк = 3^(Кк-і + Кк)/8;
ЬК(К-і) = -і/КК-і, ЬКК = і/КК-і И СК(К-і) = аКк-і /8,СКК = 3оКК-і /8 •
Учет граничных услов ий, соответствующих строкам (3),
Р(Х)х=0 = Рлев , Р(Х)х= = Рпр ;
Р (х)х=о = рлев , дР/дх |х=г = ас2^ /4;
дР/ 5х 1 х=о = ос2Члев /4, дР/дх |х=/ = ос2спр/4
требуют следующих преобразований в матрице А и векторе правой части G:
аіі = і , аКК = і,
О = (§і = Рлев, §2 = §2 - ^і^лев, §3,!, §К-і = §К-і - а(К-і)КРпр, §К = ^Лр ) , аіі = і, аКК = аКК ,
О = (§і = Рлев , §2 = §2 - а2іРлев , §3, ! 5 §К-і, §К = §К - "Чпр /4)",
аіі = аіі , аКК = аКК ,
О = (§і = §і +а-С2Члев /4,§2 = §2 -а2іРлев,§3,!,§К--К = §К-2 /4)т•
Помимо граничных условий решение уравнения (5) должно удовлетворять заданному начальному условию Р(х,ґ)\ґ = 0 = Р0(х), в частности, (4), если исходный режим является стационарным. Следовательно, учитывая, что вектор Р(0) = (Р0(хі), Р0(хК))т можно построить следующую итерационную процедуру. С использованием метода конечных объемов (МКО) строятся матрицы В и С и вектор а по ним формируются
А = В + С / Дґ и Р = СР('№> / Дґ. Затем решается уравнение АР(і) = Р и определяется вектор решения Р('Г> на первом временном слое. Далее строятся новые матрицы В и С, а по ним формируются матрица А = В + С / Дґ и вектор Р = СР('Г> / Дґ, решается уравнение АР('Т> = Р и определяется вектор решения Р(2) на втором временном слое. Продолжая этот процесс, получим решение начально-краевой задачи для уравнения (5) на всех
временных слоях ^, п = 1,[Т / ] рассматриваемого периода
времени моделирования [0, Т].
Процедура решения уравнения (6) строится аналогично, однако внешние итерации на каждом временном слое значительно упрощаются, поскольку матрица С при конечнообъемной аппроксимации имеет ненулевые элементы только на
^ п(п + 1)
главной диагонали и вектора решения Р определяются явным образом вектором решения Р('п) из предыдущего слоя. Однако для обеспечения устойчивости вычислительного процесса уже нужно соблюдать определенные соотношения между шагами по времени и пространственной координате.
3. Постановка и методы решения задачи для магистрального газопровода с промежуточными отборами газа и КС
Магистральные газопроводы с КС и промежуточными стоками (активные ветви ГТС) или без них (пассивные ветви ГТС) обычно соединяют: а) либо концевые узлы, соответствующие основным источникам или потребителям, с узлами сопряжения с закольцованными подсистемами (замкнутыми контурами), б) либо смежные узлы закольцованных подсистем. Нестационарное течение в таких ветвях газопроводах описывается уравнением типа (2), но уже с разрывными (кусочно-постоянными) коэффициентами. Начально-краевая задача, моделирующая нестационарный поток газа, требует задания граничных условий первого, второго и третьего родов в концевых и промежуточных узлах, соответствующих КС (скачок давления) и стокам со скачкообразными изменениями давления и массового расхода. Иначе говоря, задания условий баланса либо давлений Р(х + 0) - Р(х - 0) = АР(х), либо расходов
д(х + 0) - д(х - 0) = Ад(х), либо и то и другое вместе. Следовательно, решение задачи моделирования нестационарного течения газа в магистральных
237
газопроводах рассматриваемого типа возможно лишь обобщенном смысле. Например, в пространстве Соболева 1)2 зЬ2, т. е. в подпространстве один раз дифференцируемых и суммируемых с квадратами функций вместе с их первыми производными.
Сеточная аппроксимация уравнения (2) с разрывными коэффициентами, соответствующего магистральному газопроводу (активной ветви ГТС) производится следующим образом. Множество узлов сеточного разбиения должно содержать точки расположения промежуточных КС, притоков и отборов, включая ПХГ. Для конечного объема 0.к = Ок-1 и &!к, соответствующего КС или источнику (приток, отбор, ПХГ) расположенному в активном узле к, сеточная аппроксимация уравнения (2) по пространственной координате должна учитывать скачкообразное изменение давления или расхода до и после рассматриваемого узла, т. е. для активного узла с КС:
акк = 1/кк-1 +1/кк + 3°(кк-1 + кк )/8
£к = кк-1^(3Рк (х - °) + Рк-1 )/8 +
+ кк °(3Рк (х + °) + Рк+1)/ 8, а для активного узла с источником:
ак(к-1) = -1/ кк-1 + 3° к-1 кк-1 / 8,
акк = 1/кк-1 + 1/кк + 3(° к-1кк-1 +° ккк )/8,
ак(к+1) =-1/кк + 3°ккк /8,
где значения параметров ок _ 1 = ,Я|д(п)(х - 0)|/й, ок = Я|д(п)(х + 0)|/й.
4. Постановка и методы решения задачи для ГТС или ЕСГ в целом
Балансовые модели для сети газопроводов произвольной конфигурации можно создавать аналогично изложенному выше. Для этого достаточно указать способы построения конечных объемов, соответствующим узлам соединения магистральных газопроводов.
№-1
к-
>
£+1
Рис. 1. Узел к сети газопровода
В частности (см. рис. 1), для узла к сети газопровода, где соединяются три магистрали (ветви), конечный объем можно
определить соотношением Ок = Ок'_1 и Ок и ОN+1, к = 1, N.
Тогда соответствующее балансовое уравнение для конечного объема Ок вокруг узла к будет определяться следующими ненулевыми значениями элементов матрицы А и вектора правой части О
Для произвольной конфигурации ГТС вокруг каждого узла, где соединяются три и более магистральных газопровода, аналогичным образом строится конечный объем (звездообразный) и вычисляются соответствующие ненулевые элементы А и О.
Ясно, что структура матрицы А будет определяться порядком нумерации узлов сеточной аппроксимации на магистральных ветвях сети газопроводов ГТС, и возникает необходимость выбора оптимального порядка нумерации, обеспечивающего
1/Нк + 3акИк /8,
§к = 4_1(3° рк +° р-1)/8+\+ек+А+1)/8+К+1(3еЛ +с№1р+1)/8.
наибольшую эффективность балансовых математических моделей. Иначе говоря, узлы аппроксимации нужно пронумеровать так, чтобы структура матрицы А обеспечивала возможность декомпозиции вычислительной процедуры, обеспечивающей минимальные затраты времени и объема вычислений. Для достижения этой цели можно построить ориентированный граф с множеством вершин и ребер, соответствующим узлам и ветвям рассматриваемой ГТС. Для структурной декомпозиции можно воспользоваться эвристическими методами параллельных и вложенных сечений с нумерацией узлов согласно обратному алгоритму Катхилла-Макки (см. [1]). Обычно структура ГТС отличается наличием сложных подсистем с вложенными контурами, соединенных вытянутыми по протяженности (возможно, древовидными) сетями газопроводов с малой шириной вдоль полосы протяжения. Для таких систем изложенные методы обеспечивают возможность представления А в виде согласованно упорядоченной блочной матрицы с профильной стреловидной структурой
к=
0 ■ • 0
0 ф . • 0
0 0 ■ . ®
® . . ®
о о
в;
г<%> 0 ■ • 0 ф
0 ® . ■ 9
0 0 ■ . ®
® . . $
рис. 2. Структура матрицы А Таким образом, матрица А имеет следующее представление:
(Н Л , т
(9) А = * _ , Бк = diag (П*£йт) ,
^ Н У в 0
где матрицы Огк, i = 1, т, Вв и , i = 1, т - обычные или блочные тридиагональные, если соответствующая ветвь газопровода однониточная или многониточная, Н = (Н1, Н2, ..., Нт)Т и НТ = (НТ1, НТ2, ..., НТт) - разреженные коммутативные блочные матрицы связей (некоторые блоки могут быть нулевыми). При этом соответствующая структурная декомпозиция и нумерация узлов определяются соотношениями & = (&Я = & и &2 и ... и &т) и &в и
I "т)
(10) I(О) =
1(0,) I (О2)
о,..., ^1,^1 +1,..., N1 + #2,...,1 + 1 N,..., I ,...,1 + 2 N.,..., Мв + 2 N.
1(0* ) I (Ов )
I (О.) = {N.-1 +1,..., М-1 + ^, N-1 + ^ +1,..., N.-1 + N,Л + ^},
N о = 0,11 (О.)| = N.,. = 1т где I определяет нумерацию узлов, индекс в определяет нечетный узел (участок), Я - четный.
Необходимо отметить, что блочно-тридиагональное представление (9) матрицы А является наиболее эффективным при моделировании многониточных коридоров магистральных газопроводов с перемычками.
В действительности для ГТС с любой сложностью конфигурации, соответствующая матрица А может быть представлена в аналогичной (9) блочно-диагональной форме. Может измениться лишь число, расположение и тип вложенных блоков.
5. Блочное распараллеливание вычислений
В соответствии с рекомендациями [4] для решения уравнения АР = Р на и-м временном слое целесообразно использовать метод верхней релаксации по следующей схеме
т-1
ОеРЯ2= т(Ря -НРв2->) + (1 - т)
(11) ВвР{2’) = т(Рв - НТР$-") + (1 - т)ОвРв2-2\
где 1 < ю < 2 / (1+ 3) - параметр релаксации (3 - константа), верхние индексы (^ + 1) и (25) соответствуют номерам итераций по пространственной координате. Причем процедуры вычислений по первой и второй строкам (10) производятся поочередно только на нечетных и четных итерациях, соответственно.
Если матрица А имеет вид (9), то многоуровневая реализация метода реализуется согласно вычислительной схеме с тремя уровнями иерархии. На верхнем уровне иерархии по заданному начальному приближению Р(0)в и
ОкР% = ю(Ря - НР(0)в) + (1 - ю)ВйР(0)я 1
производится первая (нечетная) итерация и определяется Р(1)я. Затем производится вторая (четная) итерация
ВйР(2)в = ю(Рв - НТР%) + (1 - ю)ВвР(0)в и определяется Р(2В. Следующие парные итерации (5 := 5 + 1) производятся аналогичным способом, если выполняется условие абсолютной нормы ||Р(2? + !) - Р(2? - :)|| > е. Реализация первой итерации ОЯР(Г)Я = ю(РЯ - НР(0В) + (1 - «)ВЯР(0)Я и последующих нечетных итераций, т. е. 25 + 1, 5 > 1, состоит из аналогичных вложенных итерационных процедур для решения подзадач
( ПІ и Лґ ТУіЛ
( ) ^ н
среднего (второго) уровня
к
нт
РІ
гя
РІ
1 гг
Ср ^
гя
1 г?
І = 1, т .
Для этого с использованием формул (11) сначала определяются решения Р(г)в, . = 1, т, векторных уравнений
Вврв = Рв - НТ. Р’я, . = 1, т . Затем при этих фиксированных значениях Р(г)в, . = 1, т , определяются решения Ря, . = 1, т, подзадач нижнего (третьего) уровня иерархии, т. е. векторных уравнений ВЯР'Я = Ря - НР'в, . = 1, т . Наконец, определяется вектор Р2 + !) = ((Р:я, Р 2Я)(25 + ч (Р1в, Р2в, Р3в)(25))Т, 5 = 1, 2, ..., с использованием которого можно производить следующую итерацию (положив 5 = 5 + 1) для решения уравнения 242
АР2 + !) = Р на верхнем уровне. Таким образом, основной процесс вычислений состоит из трех вложенных идентичных процедур. Для однониточных и многониточных ветвей ГТС с обычной и блочной тридиагональными структурами для решения балансовых уравнений на любом уровне иерархии можно воспользоваться простым и очень эффективным методом прогонки.
Для повышения точности и эффективности предлагаемого подхода без значительного измельчения сетки (без увеличения размерности балансовых уравнений) необходимо воспользоваться многосеточными методами [6], основанные на схеме Р. П. Федоренко [4, 5]. В частности, распараллеленным вариантом универсальной многосеточной технологии [3] (УМТ) для одномерных по пространственной координате задач моделирования [2].
Базовый алгоритм балансового моделирования состоит из следующих этапов:
1) построение сетки & = (&я = & и &2 и ... и &т) и &в с оптимальным упорядочением узлов (10),
2) конечно-объемная (балансовая) аппроксимация краевой задачи на построенной сетке и соответствующей рассматриваемой ГТС или ЕСГ в целом,
3) построение и приведение сеточных уравнений к виду упорядочение АР = Р,
4) решение системы уравнений АР = Р с многоуровневым распараллеливанием вычислений.
6. Многосеточная схема распараллеливания вычислений
Если числа ^. = 1, т , и N'в,. = 1, т, выбраны кратными трем, можно построить иерархическую многосеточную структуру для распараллеливания вычислений, представленную на рис. 2. Множество сеточных точек балансовой аппроксимации любой ГТС представляется как объединение множества узлов и
границ конечных объемов линейных участков, т. е. нулевой уровень грубости [2]
0(0,1) = О у(0,1) и0ґ(0,1);
Оу (0,1) = {х1 : х1 еП, \ = хк+1 - х1, к = 0, N;
Ог (0,1) = {хк : хк = 0,5(хк + хк+1), к = 1^, где N - количество точек разбиения.
Исходную сетку &(0, 1) представим как объединение трех более грубых и не пересекающихся сеток первого уровня, т. е.
3
П(0,1) = и О(1,«), П(1,а) пП(1,Ь) = 0, аФ Ь
а=1
Далее рекуррентным способом каждая из сеток П(1,а),а = 1,3, рассматривается как исходная для сеток П(2,а),а = 1,!,32, а полученные девять еще более грубых сеток образуют второй уровень т. д. Построение грубых сеток заключается в удалении двух точек из & и &, как показано на рис. 3. Построенная таким образом иерархическая многосеточная структура с тремя уровнями грубых сеток будет иметь вид, представленный на рис. 4.
0(0,1)
0(1,1)
0(1,2)
0(1,3)
X* у г
* ч / *
* I • I * I * I » I « I *-1 - 1*1* I * ^
I
‘1‘
..
А
I '
і
Ц!
.
А..
11
I
і.
і
.....\
0-^ уровень
£
)И
Рис. 3. Построение грубых сеток
Нулевой уровень
Рис. 4. Многосеточная структура
Многосеточная схема иерархических вычислений может быть организована следующим образом. Исходные системы уравнений типа (7) или (8), соответствующие рассматриваемой ГТС или ЕСГ, заменой переменных Р = Р + и преобразуются к
(7 ' ) (— С+в)(Р+и)(п+1) = — СР(п) ^(—С+ви(п+1) = — СР(п) -(— С+в)Р("+1),
А/ А? А? А/ А/
(8 ' ) (— С +Б)(Р+и)(и+1) = — СР(и) ^(— С +5)и(и+1) = — СР-(— С +Б)Р(п+1\
А/ А А/ А/ А/
С учетом граничных условий типа (3), преобразованных аналогичным образом, системы (7 ) или (8 ) можно преобразовать и получить систему уравнений нулевого уровня
(11) А 0и 0 = J 0, J 0 = (^ 1, ^ N ) , и0 = (и1, !, UN ) ,
где матрица А0 имеет характерную блочно-стреловидную структуру, приведенную выше. Системы уравнений первого, второго и т.д. уровней строятся согласно схемам, представленным на рис. 3 и 4. Для трехуровневой многосеточной структуры уравнения, соответствующие узлам сетки третьего уровня, могут быть построены согласно схеме на рис. 5, т. е.
4
ч
* I
*
4
*10
'{1} Ч:
О-й уровень (сам ія мелкая сетка)
„V
4
{2}
■■И"
■14
{2}
3-й уровень (самая грубая сетка)
Рис. 5. Построение конечных объемов на грубых сетках где хУ{2} определяет соответствующую узлу {2} строку матрицы
А 3. ^ з
Интеграл J = {2} в правой части уравнения, соответствующего конечному объему [Х{1}, Х{2}] ° [Х5, х^], на грубой сетке третьего уровня должен вычисляться как сумма интегралов по девяти конечным объемам на самой мелкой сетке нулевого уровня, поскольку [./{1}, Х/{2}] °
° [/5, /м] = [/5, /6] и [/;, /7] и ... и [Х^, /13] и [Х^з, /м].
Строка матрицы А3 , соответствующая конечному объему [х^{!}, х/{2}], вычисляется непосредственно на грубой сетке третьего уровня обычным способом. Аналогичным образом могут быть построены системы уравнений для всех уровней
(12) А,и, = Зь, ^ = (J1І,..., JLNl)т,
и 1 = (и1,...,п11)т, I = 1,.,I.
где 1+ - номер уровня с самой грубой сеткой. Многосеточные итерации производятся согласно схеме, изображенной на рис. 6. Компоненты вектора и1, вычисленные на грубых сетках, добавляются к соответствующим компонентам вектора приближенного решения Р, полученного на предыдущем шаге итерационной процедуры.
\
2
3
4
5
Рис. 6. Схема распараллеливания на грубых сетках
Если временем пересылки данных между основным процессором и процессорными модулями можно пренебречь, то асимптотические значения ускорения Бм = 7(1) / Т(М), где 7(1) и Т(М) - затраты времени при использовании одного процессора и М процессоров, и эффективности Ем = Бм / М, согласно [3] определяются соотношениями
Бм = М - (М - 1)Т0/а , Ем = 1 - (М - 1)Т0/М,
/ 1=1 / 1=1
где Т0 и Т1 - затраты времени на нулевом и 1-м уровнях. Ускорение и эффективность могут быть улучшены, если в процессе решения задачи результаты вычислений на грубых сетках первого уровня принять как базовые, которые принято называть динамическими [3], а пилообразный динамический цикл выполнять по схеме, представленной на рис. 7.
Для конкретных ГТС или ЕСГ в целом пилообразная многосеточная схема вычислений с динамическим циклом используется для каждого уровня 1 < I < 1+, начиная с самых грубых сеток на уровне 1+, т. е., согласно рис. 7, если 1+ = 5. Другими словами, решаются соответствующие системы уравнений вида
(12) для любой грубой сетки Щ(Ь, 1), Щ(Ь, 2) и Щ(Ь, 3), если 1 < Ь < Ь+.
Таким образом, для каждого временного слоя балансовые модели ГТС и ЕСГ в целом определяются решениями уравнений вида (12) , а сам процесс вычислений состоит из Ь+ + 1 уровней иерархии многосеточной схемы. На любом уровне многосеточной схемы для всех грубой сеток этого уровня иерархии процесс вычислений производится с использованием циклической схемы (10) согласно структуре разбиения матрицы А0, т. е.
БД, / = 1, т , Ов и БВ, / = 1, т . Следовательно, каждый такой фрагмент внутренней процедуры вычислений, в свою очередь, состоит из трех уровней иерархии блочных циклических итераций вида (10).
Сходимость многоуровневых вычислений при балансовых способах моделирования нестационарных потоков газа в сетях газопроводов любой сложности определяется сходимостью решений последовательности линейных систем уравнений по пространственным координатам при фиксированных значениях
о = 2Л|д1п)| / й, п = 1, 2, ... на каждом временном слое. Многократные вычислительные эксперименты на модельных задачах для сложных сетей газопроводов подтверждают сходимость этих последовательностей к обобщенному решению нелинейной системы. Однако строгое доказательство сходимости этой последовательности требует дополнительных исследований.
7. Заключение
Основными достоинствами балансовых моделей распределения давления и потоков газа в сетях газопроводов произвольной конфигурации является их предельная простота и универсальность. Эти особенности в наибольшей степени проявляются при учете структурных особенностей сетей газопроводов любой сложности и соответствующих сеточных аппроксимаций распределения давления и потоков газа в них с точностью, практически определяемой погрешностью исходных данных. Простота и универсальность таких моделей проявляется также при учете влияния управляющих воздействий в промежуточных узлах, где расположены КС и другие объекты (источники, ПХГ, буферные и др. потребители) с регулируемыми подачами и отборами газа.
Использование внешних иерархических многосеточных итераций с пилообразным динамическим циклом на всех грубых сетках каждого уровня Ь : 1 < Ь < Ь+, в сочетании циклическими внутренними итерациями (10), обеспечивает почти не улучшаемые показатели эффективности предлагаемых вычислительных методов моделирования.
Балансовые модели нестационарных режимов течения газа в сетях газопроводов произвольной конфигурации предоставляют возможность решения основных задач диспетчерского управления ГТС и ЕСГ в целом в наиболее общей и актуальной постановке.
Литература
1. ДЖОРДЖ А., ЛЮ ДЖ. Численное решение больших разреженных систем уравнений М.: Мир, 1984.
2. МАРТЫНЕНКО С.И. Универсальная многосеточная технология для численного решения дифференциальных уравнений в частных производных на структурированных сетках // Вычислительные методы и программирование. 2000. 1, № 1. 83 - 102.
3. МАРТЫНЕНКО С.И. Распараллеливание универсальной многосеточной технологии // Вычислительные методы и программирование. 2003. 4, № 1. 49 - 54.
4. ФЕДОРЕНКО Р.П. Релаксационный метод решения разностных эллиптических уравнений // Журнал вычислительной математики и математической физики. 1961. 1, № 5. 992 -927.
5. ФЕДОРЕНКО Р.П. Итерационное решение разностных эллиптических уравнений // Успехи математических наук. 1964. 28, вып. 2. С. 121 - 182.
6. HACKBUSCH W. Multigrid Methods and Applications. Berlin, 1985.
MULTIGRID BALANCE MODELS OF UNSTEADY FLOWS IN COMPLEX GAS TRANSPORTATION SYSTEMS
Atlas Akhmetzyanov, Institute of Control Sciences, RAS, Moscow, Cand. Eng. Sc. (Moscow, Profsoyuznaya st., 65, (495) 334-90-30, awa@ipu.ru).
Sergey Spiridonov, Institute of Control Sciences, RAS, Moscow (spirid@ipu.ru).
Anton Salnikov, Institute of Control Sciences, RAS, Moscow (salnikov@ipu.ru).
Abstract: Multigrid variants are suggested for hierarchical flow distribution models to optimize unsteady conditions of gas transport at various levels of dispatcher control of the Unified Gas Supply System of Russia. The multigrid scheme based on the relaxation scheme by R. P. Fedorenko is used to increase calculative accuracy of multilevel methods for each level of the hierarchy of calculations. Generality of the approach is determined by the ease of construction of a multigrid structure and transition operators that provide high-precision discretization on coarse grids without preliminary smoothing and interpolation inherent in classical multigrid methods. The universal multigrid technique in the finite-volume discretization is suitable for objects described by equations with discontinuous and nonlinear coefficients and right members.
Keywords: multigrid method, balance model, unsteady flow, gas transportation.
Статья представлена к публикации членом редакционной коллегии М. В. Губко
