Гусева А.Б., Кушнаренко В.Г. ДИСКРЕТНАЯ SIMULATION-МОДЕЛЬ СЕТЕВОГО ДИНАМИЧЕСКОГО ОБЪЕКТА С РАСПРЕДЕЛЕННЫМИ ПАРАМЕТРАМИ НА ОСНОВЕ БЛОЧНОГО ЧИСЛЕННОГО МЕТОДА

ДИСКРЕТНАЯ SIMULATION-МОДЕЛЬ СЕТЕВОГО ДИНАМИЧЕСКОГО ОБЪЕКТА С РАСПРЕДЕЛЕННЫМИ ПАРАМЕТРАМИ НА ОСНОВЕ БЛОЧНОГО ЧИСЛЕННОГО МЕТОДА

А.Б. Гусева (аспирант)
В.Г. Кушнаренко (магистрант)
Донецкий национальный технический университет
guseva.annab@gmail.com
kushnar_vova@mail.ru

Рассматривается формальная постановка задачи моделирования сетевого динамического объекта с распределенными параметрами (СДОРП) и преобразование СДОРП-модели на основе блочного численного метода (БЧМ) к форме дискретной Simulation-модели (ДSМ), являющейся базой для построения параллельного решателя уравнений. Предлагаются расчетные формулы ДSМ, основанные на одношаговом четырехточечном БЧМ.

Математическое моделирование, сетевой динамический объект, блочный численный метод, дискретная Simulation-модель.

Введение

Сложные сетевые динамические объекты с распределенными параметрами (СДОРП) различных предметных областей требуют при их исследовании высокопроизводительных средств и методов решения многомерных систем уравнений, а также компьютерной поддержки всех этапов разработки и применения средств моделирования [1, 2]. Модель СДОРП – это формальное представление топологии и система дифференциальных уравнений в частных производных, описывающих динамические процессы в ветвях и узлах объекта. Simulation-моделью (SM) принято называть модель, преобразованную к виду, удобному для применения численных методов. Дискретную Simulation-модель (DSM) получают посредством аппроксимации SM-уравнений, граничных и начальных условий расчетными формулами выбранного численного метода. Аппаратно-программная реализация DSM представляет собой решатель системы уравнений модели. СДОРП-симулятор реализует все этапы работ по моделированию – от представления топологии объекта до проведения модельных экспериментов. Перспективным направлением является разработка СДОРП-симуляторов в составе распределенных параллельных моделирующих сред (РПМС) общего и проблемно-ориентированного назначения [2]. При этом наряду с известными численными методами, применяемыми в последовательных языках моделирования и реализуемыми в РПМС, требуется построение параллельных решателей на основе новых блочных численных методов (БЧМ), в которых для каждого блока, состоящего из k точек, новые k значений функции рассчитываются одновременно [3]. Рассмотрим ключевой этап в разработке параллельных БЧМ-решателей СДОРП-симуляторов – преобразование СДОРП-модели к ее DSM-форме.

Формальное описание (модель) и постановка задачи параллельного моделирования СДОРП

Топология аэродинамического сетевого динамического объекта характеризуется графом G (m, n), который кодируется таблицей 1, имеющей m строк и s+5 столбцов. Здесь QJ – поток воздуха в ветви j; AKk и EKk соответствуют начальному и конечному узлам j-ветви; kÎ (1,2,..., n), jÎ (1,2,..., m); PAR(PJ1, PJ2, …, PJs) – множество s параметров PJ каждой ветви; AEJ – активный элемент в j-й ветви; VECOMJ – вербальный комментарий, объясняющий технологическое назначение j-ветви.

Таблица 1 - Исходное кодирование графа G (m, n)

AKk

EKk

PAR(PJ1, PJ2, …, PJs)

AEJ

VECOMJ

Для j-й ветви, не имеющей утечек воздуха, динамика расхода и давления описываются системой уравнений [4]

, (1)

в которой Pj, Qj – давление и расход воздуха вдоль координаты x, отсчитываемой от начального AKk до конечного EKk узлов; rj – удельное аэродинамическое сопротивление; Fj – площадь поперечного сечения ветви; r – плотность воздуха; a – скорость распространения звука в воздухе; rj(ξp,t) – регулируемое сопротивление; ξp – координата регулирующего органа.

Граничные условия для системы уравнений (1) – это функции давления в начальном РAKk и конечном PEKk узлах ветви i. По типу граничных условий ветви сети делятся на три вида:

Ветви, инцидентные внутренним узлам сети, в которых давление рассчитывается в процессе решения сетевой системы уравнений в соответствии с узловыми динамическими условиями

, (2)
где Pwj – давление в узле WJ; Qwj – суммарный поток воздуха в узле WJ; Fwj – площадь поперечного сечения j-й ветви в узле;

Ветви, инцидентные узлу подключения вентилятора; в этом узле давление задается как характеристика вентилятора

Pwj = PAEj(Qj); (3)

Ветви, инцидентные узлу выхода в атмосферу; здесь принимается постоянное давление атмосферы

Pwj = PAТМ = const. (4)

Для всего графа имеем n = n1 + n2 + n3 граничных условий, где n1, n2, n3 – количества узлов с условиями (2), (3), (4) соответственно.

Начальные условия – заданные значения переменных вдоль координаты ξ:

Pj(ξ, 0), Qj(ξ, 0) (j = 1,2,...,m) (5)

Задачу параллельного моделирования СДОРП сформулируем следующим образом: для объекта, кодируемого таблицей 1 графа G (m, n) и описываемого m системами уравнений (1) с n граничными (2), (3), (4) и 2m начальными условиями (5), разработать и имплементировать параллельные алгоритмические и аппаратно-программные средства, которые адекватно воспроизводят динамические процессы Pj(ξ, t), Qj(ξ, t) (j = 1,2,...,m) при управлении расходами воздуха с помощью изменения граничных условий (3) и аэродинамических сопротивлений rj(ξp,t), а также действии возмущений производственного характера.

Дискретизированные по пространственной координате Simulation-модели j-й ветви и всего СДОРП

В результате аппроксимации системы уравнений (1) по методу прямых с шагом Δξ по пространственной координате получим следующие уравнения для k-ого элемента j-ой ветви:

(6)

Внутренние граничные условия типа (2) аппроксимируются уравнением:

(7)

Здесь Pwi = PjMj+1 – давление в конечном узле последнего элемента j-ветви с расходом Qjk = QjMj, который направлен в wi-узел; jwiЄj– номер ветви из множества j=1, 2,..., m, которая инцидентна узлу wi; при этом jwi1 является первым элементом j-ветви с wi как начальным узлом (выток), а jwiMj – конечным элементом с wi как конечным узлом (приток). Каждая ветвь разделяется при аппроксимации на Mj элементов Qj1,...,QjMj. При этом значения давлений нумеруются как Pj1, Pj2,...., PjMj+1. Следует отметить, что j-ветвь имеет начальный узел wi с давлением Pwi=Pj1 и конечный узел wi+a (a = const) с давлением Pwi+a = PjMj+1.

Для аэродинамического объекта с распределенными параметрами каждая ветвь представляется в виде двух векторов Qj, Pj (j = 1…m):

Qj = (Qj1,Qj2,…,QjMj)T (8)

– поток воздуха в j-ой ветви;

Pj = (Pj1,Pj2,…,PjMj+1)T (9)

– давление в j-ой ветви.

При этом количество элементов в векторах

Mj = lj / Δξ (10)

вычисляется по длинам ветвей lj и единому для СДОРП шагу Δξ.

Дискретизированная по пространственной координате модель j-ветви содержит Mj систем уравнений (6) при k=1,2,...,Mj. Для вычисления компонент векторов (8), (9) приведем систему (6) к форме пространственно-дискретизированной Simulation-модели (ПДSM)

(11)

с коэффициентами αj, βj, βrj, gj, зависящими от параметров j-ветви.

При разработке ПДSM для всего СДОРП необходимо представить mсистем уравнений типа (11) для всех ветвей при j=1,2,...,m:

k=1,2,...,M1

::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: (12)

k=1,2,...,Mm

В соответствии с (2) сформулируем n1граничных условий для внутренних узлов СДОРП (wi=1,2,...,n1):

:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: (13)

Характеристики активных элементов (3) дают n2граничных условий вида

Pwi = P1AEJ(Qj)

:::::::::::::::::::::::: (14)

Pwi = Pn2AEJ(Qj).

При этом wi=n1+1,n1+2,..., n1+n2 , если нумерация узлов продолжает номера из (13).

Для n3 узлов подключения к атмосфере имеем по (4):

Pwi =P(n1+n2)+1= P1AТМ = const,

:::::::::::::::::::::::::::::::::::::::::::: (15)

Pwi=Pn= Pn3AТМ = const.

При указанной нумерации имеем в (15) wi=(n1+n2)+1,...,(n1+n2)+(n3-1), n.
Система (12) вместе с граничными условиями (13), (14), (15) образует пространственно-дискретизированную СДОРП-Simulation-модель.

Построение дискретной модели СДОРП на основе блочного численного метода

Для задачи Коши

(16)

временная ось, состоящая из М дискретных значений, разбивается на N блоков с шагом τ, каждый из которых состоит из k точек (рис 1). Каждая нулевая точка (tn,0) является начальной для расчета следующих точек блока, но при этом не входит в него.

Рисунок 1 – Схема разбиения временной оси на блоки

Существует два типа блочных методов: одношаговые и многошаговые многоточечные [2, 3]. Для одношаговых блочных методов характерно определение начальной точки блока, исходя из одного рассчитанного значения последней точки предыдущего блока. Если для расчета нулевой точки n-ого блока принимают участие две и более точки предыдущего блока, речь идет о многошаговом блочном методе.

В общем виде уравнение m-шагового k-точечного блочного метода записывается следующим образом [3]:

(17)

Здесь:

un,i(iτ)– приближенные решения в точке i блока n в момент времени iτ; n=1,2,...,N; i = 1,2,...,K.

un,0 – приближенное решение в нулевой точке блока n; для первого блока данная величина будет являться начальным условием задачи Коши, для последующих блоков нулевая точка каждого блока будет совпадать с последней точкой предыдущего, т.е. u2,0 = u1,K,..., uN,0 = uN-1,K;

Fn,j = f (tn+jτ; un,j) – значение правой части в точке j блока n, где j = 1,2,...,m;

Fn-1,j = f (tn-1+jτ; un-1,j) – значение правой части в точке j блока n;

ai,j, bi,j – коэффициенты, предварительно рассчитываемые с помощью построения интерполяционного многочлена Лагранжа.

Преобразуем пространственно-дискретизированную СДОРП-Simulation-модель (12), (13), (14), (15) в ДSМ-форму на основе одношагового (m=1) многоточечного блочного численного метода (рис.2)

Рисунок 2 – Шаблон одношаговой k-точечной разностной схемы.

Формулы для параллельного двухточечного блочного метода можно получить из (17), предварительно рассчитав коэффициенты ai,j и bi с помощью построения интерполяционного многочлена Лагранжа [2, 3].

При k=4 расчетные формулы блочного метода будут иметь такой вид [3]:

(18)

Возвращаясь к системе уравнений (12), изменим индексацию элементов ветвей с целью корректного представления дискретной модели:

(19)

Здесь v – номер элемента j-ой ветви (v=1,2,...Mj). Следует заметить, что слагаемое βrjQjv|Qjv| соответствует регулируемым сопротивлениямrj(ξp,t), которые устанавливаются в одной ξp-точке заданного v-элемента j-ветви.

Предлагается следующая методика разработки ДSM для j-ветви:

Выбираются параметры j-й ветви из табл. 1.
Задается Tsim – длительность исследуемых динамических процессов, сек.
Определяется шаг Δξ, единый для всех ветвей СДОРП; так, для шахтных вентиляционных сетей представляет интерес диапазон шага 5.0 ≤ Δξ ≤ 50.0м.
Вычисляются шаг по времени τ=Δξ/a, количества шагов M = Tsim / τи блоков N = M/k(мы рассматриваем метод сk=4).
Для наглядного представления будущих вычислений разбивается ось времени в диапазоне 0 ≤ t ≤ Tsimпо шаблону рис.2 для k=4 на блоки с номерами n = 1, 2,..., r-1, r, r+1,..., N-1, N. Этот шаблон следует сделать для обоих уравнений системы (19), чтобы представить возможные взаимосвязи и необходимый обмен данными между Q- и P-процессами.
Согласовать обозначения переменных, коэффициентов и индексов уравнений ДSM (12), (13), (14), (15) с обозначениями в (17), (18) в виде Qjvn,i = Un,i, Pjv+1, n,i = Un,i, а в правых частях Fn,0, Fn,1 – Fn,4 сделать пометки о принадлежности их к Qj- или Pj-уравнениям.
Записываются следующие системы БЧМ-ДSM-уравнений для вычисления Qj- и Pj-значений в каждой точке n-го блока:

(20)

Значения Qvn,0 и Pv+1,n,0 в нулевой точке блока n приравниваются значениям в последней точке предыдущего блока для n >1 (tn-1,k = tn,0, рис.2). Правая часть Fn,k для нахождения переменных Qjv из (19) вычисляется по следующим формулам:

Функция Fn,k для вычисления Pjv+1 из (19) рассчитывается как:

(Fv,n,0)Pjv+1 = gj(Qjv,n,0 – Qjv+1,n,0)
(Fv,n,1)Pjv+1 = gj(Qjv, n,1 – Qjv+1,n,1)
(Fv,n,2)Pjv+1 = gj(Qjv, n,2 – Qjv+1,n,2) (22)
(Fv,n,3)Pjv+1 = gj(Qjv, n,3 – Qjv+1,n,3)
(Fv,n,4)Pjv+1 = gj(Qjv, n,4 – Qjv+1,n,4)

СДОРП-ДSM получим, распространив предложенную методику на все ветви и узлы:

Для каждой СДОРП-ветви в соответствии с (12) при j=1,2,...,m составляется БЧМ-ДSM (20) – (22).
Внутренние граничные условия (13) представляются по БЧМ-формулам (18) с учетом обозначений, используемых в (20), (21), (22).
Внешние граничные условия (14), (15) придаются соответствующим инцидентным ветвям СДОРП.

Ввиду большой размерности уравнений БЧМ-ДSM сетевого объекта их формирование следует реализовать как функцию генератора уравнений РПМС.

Выводы

Для параллельного моделирования сложных сетевых динамических объектов с распределенными параметрами необходима алгоритмическая и программная поддержка построения решателей уравнений. В статье предложена и проиллюстрирована методика преобразования исходной модели СДОРП в дискретную форму на основе блочного четырехточечного численного метода. Этапы преобразований служат базой для реализации генератора уравнений виртуальной параллельной модели и MIMD-симулятора в составе распределенной параллельной моделирующей среды.

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