Колодкин В.М., Мурин А.В., Петров А.К., Горский В.Г. Количественная оценка риска химических аварий

Подождите немного. Документ загружается.

5.4. Подсистема «Атмосфера» 121

Разностная сетка и метод расщепления. В настоящей работе при числен-

ном решении уравнение (4.4.4) аппроксимировалось неявной консервативной раз-

ностной схемой, заданной на неравномерной прямоугольной разностной сетке. При

построении разностной схемы использовался метод расщепления по физическим

процессам и пространственным направлениям в сочетании с методом контрольно-

го объема [176, 177]. Структура разностной сетки проиллюстрирована на рис. 5.6.

Большинство физических параметров относятся к центрам соответствующих ячеек

сетки. Коэффициенты обмена и потоки относятся к граням ячеек. Граничные усло-

вия выставляются на соответствующих гранях ячеек. Такой подход обеспечивает

консервативность схемы, что особенно важно при моделировании распростране-

ния аварийных воздействий. Использование прямоугольных контрольных объемов

способствует простоте программной реализации, а неравномерность (сгущение в

областях больших градиентов) сохраняет приемлемую точность.

Для перехода с дискретного временного слоя n на следующий временной слой

используется следующая схема расщепления:

n+1

= (T

· f

, t = nτ, (5.4.1)

где f = c

, θ, w — переносимая физическая величина, T

— разностные операто-

ры адвекции и диффузии, соответственно, S — источниковый оператор (локальные

изменения концентрации за счет химических реакций, эмиссии и т.д.), τ — вре-

менной шаг интегрирования модели.

Дискретизация оператора адвекции-диффузии. Введем понятие суммарно-

го потока (адвективного и диффузионного) через грань контрольного объема в

направлении оси x

J = u

f − K

. (5.4.2)

Рассмотрим контрольный объем i, грани контрольного объема обозначим индекса-

ми −1/2 и +1/2 для левой и правой грани, соответственно (см. рис. 5.6).

После интегрирования по контрольному объему уравнения диффузии матери-

альный баланс для величины f в одномерном случае запишется в виде:

n+1

− f

= J

n+1

−1/2

− J

n+1

+1/2

, (5.4.3)

где h

— длина грани i-го контрольного объема. Отсюда дискретный оператор

адвекции-диффузии можно представить в виде:

(T D)

· f

= f



n+1

−1/2

− J

n+1

+1/2



. (5.4.4)

Исследование одномерной стационарной задачи показало, что точное решение

обеспечивает экспоненциальная схема [169]:

n+1

+1/2

= (u

)

+1/2

n+1

− f

n+1

i+1

exp(P e

+1/2

) − 1

, (5.4.5)

где P e

+1/2

= (u

)

+1/2

(δx

)

+1/2

/(K

)

+1/2

– число Пекле, (δx

)

+1/2

— расстояние

между центрами ячеек i и i + 1.

122 Компьютерная система прогноза последствий аварий

Расчет экспонент, входящих в выражение (5.4.5), является весьма трудоем-

кой процедурой, поэтому в программной реализации имеет смысл использовать

схему со степенным законом [169], обеспечивающую приемлемую аппроксимацию

экспоненциальной схемы:

P e

exp(P e) −1











−P e, P e < −10,

(1 + 0.1P e)

− P e, −10 ≤ P e < 0,

(1 − 0.1P e)

, 0 ≤ P e ≤ 10,

0, P e > 10.

(5.4.6)

Представленная неявная схема абсолютно устойчива и имеет второй порядок

точности по пространственным координатам.

5.4.4. Алгоритм параллельных вычислений. Естественным подходом для

использования параллельных вычислений при решении упрощенного уравнения

турбулентной диффузии (без учета поперечной диффузии) является декомпозиция

расчетной области на вертикальные области (колонки), т. е. разбиение трехмерной

сетки вдоль x и y осей и введение соответствия между физическими подобластями

и узлами двумерной логической сети процессоров [93]. Такое разбиение удобно

использовать при использовании явных численных транспортных схем [178–180].

Специфика задач по оценке краткосрочных (порядка нескольких часов) по-

следствий аварий такова, что линейный масштаб переноса вдоль направления сред-

него ветра значительно больше, чем поперечный и вертикальный масштабы. В

связи с последним обстоятельством, возможно ввести одномерное разбиение вдоль

направления среднего (за рассматриваемый период времени) ветра [11]. Отметим,

что алгоритм параллельных вычислений разрабатывался для сетевой вычислитель-

ной системы, т. е. для системы компьютеров, связанных относительно малопроиз-

водительной локальной сетью. Для вычислительных систем другой архитектуры

более эффективным может быть использование других алгоритмов.

Для распараллеливания вычислений в выражение для потока на грани ячей-

ки (5.4.5) вводится «степень явности», определенная следующим образом:

n+1

+1/2

= (u

)

+1/2

n+1

− f

i+1

exp(P e

+1/2

) − 1

. (5.4.7)

n+1

−1/2

= (u

)

−1/2

n+1

i−1

n+1

i−1

− f

exp(P e

−1/2

) − 1

Такая постановка обеспечивает, во-первых, возможность вычислять значения кон-

центраций в слое i + 1 явным образом, не прибегая к методу прогонки. Во-вторых,

схема остается абсолютно устойчивой и консервативной.

Введем следующие обозначения: P — число подчиненных процессоров (ПП),

УП — управляющий процессор, n

— число узлов разностной сетки вдоль оси

, приходящихся на ПП

, а i

∈ [0, n

− 1] — индекс, нумерующий узлы в ПП

Тогда параллельный алгоритм может быть представлен следующим образом (для

каждой итерации по времени):

1. Каждый ПП

решает систему f

n+1/6

ijk

= (S)

· f

ijk

в своей подобласти;

5.4. Подсистема «Атмосфера» 123

...

PSfrag replacements

ПП

УП

запрос

ответ

n+4/6

= 0

n+1

= N

− 1

Рис. 5.7. Схема потоков данных при параллельном решении уравнения турбулент-

ной диффузии

2. Каждый ПП

решает систему алгебраических уравнений:

n+2/6

ijk

= (T

)

n+1/6

· f

n+1/6

ijk

в своей подобласти (трехточечная прогонка);

3. Каждый ПП

решает систему f

n+3/6

ijk

= (S)

n+2/6

·f

n+2/6

ijk

в своей подобласти;

4. Каждый ПП

решает систему алгебраических уравнений

n+4/6

ijk

= (T

)

n+3/6

· f

n+3/6

ijk

в своей подобласти (трехточечная прогонка);

5. Если p 6= 0, передать в ПП

p−1

значения f

n+4/6

для слоя i

= 0;

6. Если p 6= P − 1, принять от ПП

p+1

значения f

n+4/6

p+1

для слоя i

p+1

= 0;

7. Если p 6= 0, принять от ПП

p−1

значения f

n+1

p−1

для слоя i

p−1

= N

p−1

− 1;

8. Каждый ПП

решает систему f

n+5/6

ijk

= (S)

n+4/6

·f

n+4/6

ijk

в своей подобласти;

9. В каждом ПП

решить систему f

n+1

ijk

= (T

)

n+5/6

· f

n+5/6

ijk

(потоки на

гранях ячеек определяются выражением (5.4.7), решение в каждом слое i

получается явно из решения в слое i − 1);

10. Если p 6= P −1, передать в ПП

p+1

значения f

n+1

для слоя i

= N

− 1.

Последовательность шагов 5-10 представляет собой схему конвейерного типа:

решение в каждом слое i получается из решения в слое i − 1. Схема потоков

данных, возникающих при параллельном решении, показана на рис. 5.7. В про-

цессе решения осуществляется взаимодействие только между подчиненными про-

цессорами. Потоки данных от управляющего процессора осуществляют передачу

начальных данных и сбор результатов расчета.

Отметим, что выполнение шага 5 может быть начато сразу после того, как

получены значения на временном слое n + 4/6 в слое i

= 0, а шаг 6 осуществлен

124 Компьютерная система прогноза последствий аварий

путем неблокирующего приема данных, одновременно с проведением вычислений

в слоях i

∈ [0, N

− 2], что значительно улучшает эффективность алгоритма.

Оценка параллельной эффективности. Мерой эффективности любого па-

раллельного алгоритма является коэффициент повышения быстродействия σ, опре-

деляемый как [181]

σ =

, (5.4.8)

где τ

— время решения задачи на i процессорах и коэффициент эффективности

µ = 100σ/P. (5.4.9)

Пусть размеры разностной сетки составляют n

× n

ячеек, пропускная

способность передачи данных в локальной сети v (бит/с), l (бит) — размер машин-

ного слова для хранения значений f, γ — время (с), затрачиваемое на вычисления,

приходящиеся на один контрольной объем. Будем рассматривать топологию ЛВС

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

(сетевой концентратор).

Тогда время, затрачиваемое на вычисления последовательной программой на

одной итерации, определяется выражением:

= n

γ. (5.4.10)

Время, необходимое для расчета одного временного шага параллельной програм-

мой, можно представить в виде:

= τ

+ τ

, (5.4.11)

где τ

— время синхронизации и обмена данными, τ

— среднее время, затрачи-

ваемое на вычисления на одном процессоре.

В случае равномерного разбиения числа узлов сетки, приходящихся на один

процессор, τ

будет определятся следующим выражением:

= n

γ/P. (5.4.12)

В свою очередь, время синхронизации и передачи данных можно оценить как

время, необходимое на передачу всего объема данных, используемых для обнов-

ления значений f на границах подобластей, через локальную сеть:

= 2ln

(P − 1)/v. (5.4.13)

Отсюда, коэффициент сокращения времени выполнения равен:

ρ =

2l(P − 1)

γvn

, (5.4.14)

а коэффициент повышения быстродействия —

σ =

1 + 2l(P − 1)P/(γvn

)

. (5.4.15)

Анализ выражения (5.4.14) показывает, что на увеличение параллельной эф-

фективности (сокращение времени счета) влияют такие факторы, как увеличение

5.4. Подсистема «Атмосфера» 125

0 5 10 15 20 25 30 35 40 45 50

PSfrag replacements

σ(P )

Рис. 5.8. Влияние топологии ЛВС (v = 10

) и размерности задачи на масштабиру-

емость алгоритма: 1 — идеальная масштабируемость; 2 — топология звезды, n

50; 3 — топология звезды, n

= 150; 4 — топология кольца, n

= 50; 5 — топология

кольца, n

= 150

пропускной способности сети и увеличение разрешения сетки вдоль направления

преобладающего потока.

Численные эксперименты по определению эффективности предложенного алго-

ритма были осуществлены на гетерогенном кластере, состоящим из персональных

компьютеров под управлением ОС Linux, соединенных локальной сетью Fast Eth-

ernet с пропускной способностью 100 Мбит/с (v = 10

). В качестве узлов кластера

использовались серийные одно- и двухпроцессорные персональные компьютеры с

процессорами Intel Pentium III 500-800 MHz, AMD K7 500 MHz. Для хранения

дискретных значений f

ijk

использовали тип float языка C++ (l = 64).

Решалась тестовая задача определения концентрации от мгновенного точечно-

го источника на разностной сетке 150 × 50 × 50 узлов. Число итераций — 360.

Отметим, что при измерении времени выполнения в полное время решения за-

дачи входит время на инициализацию начальных данных. Выполненные оценки

показали, что параметр γ может быть оценен значением ≈ 2 · 10

−6

Влияние топологии локальной вычислительной сети на производительность

алгоритма согласно теоретическим оценкам показана на рис. 5.8. Отметим, что на

производительность негативное влияние оказывает звездообразная топология сети.

В случае топологии в виде одномерной цепочки (или кольца), эффективность будет

существенна повышена в силу исчезновения вклада P − 1 в выражениях (5.4.14)

и (5.4.15).

126 Компьютерная система прогноза последствий аварий

0 2 4 6 8 10 12 14 16

PSfrag replacements

σ(P )

Рис. 5.9. Влияние размерности задачи и производительности ЛВС на масштабиру-

емость алгоритма: 1 — идеальная масштабируемость; 2 — v = 10

, n

= 50, расчет;

3 — v = 10

, n

= 150, расчет; 4 — v = 10

, n

= 50, расчет; 5 — v = 10

, n

= 150,

расчет; 5 — v = 10

, n

= 150, измерение

В табл. 5.1 показаны время выполнения, коэффициент повышения быстродей-

ствия и коэффициент эффективности алгоритма в зависимости от числа процес-

соров в вычислительном кластере (звездообразная топология ЛВС). Отметим, что

при увеличении числа процессоров с 1 до 8 время расчета понижается более чем в

7 раз, а параллельная эффективность порядка 90%, причем результаты получены

при достаточно неоптимальной топологии сети.

На рис. 5.9 графически представлена зависимость коэффициента ускорения от

числа процессоров, размерности разностной сетки и пропускной способности ЛВС

(звездообразная топология). Некоторый разброс измеренных значений σ от рас-

считанных объясняется, во-первых, тем обстоятельством, что задействованные в

расчете процессоры несколько отличались по производительности. Во-вторых, при

выводе зависимости (5.4.15) не учитывали возможность временного перекрытия

между расчетом и обменом данными, реализованной в программе.

Таблица 5.1. Зависимость характеристик эффективности алгоритма от числа

процессоров в вычислительном кластере.

P 1 2 3 4 5 6 7 8

, с 236.0 124.0 88.0 67.0 55.6 49.0 37.2 32.2

σ — 1.9 2.7 3.5 4.3 4.8 6.3 7.3

µ, % — 95 90 88 86 80 90 91

5.4. Подсистема «Атмосфера» 127

5.4.5. Тестирование численного алгоритма. Задача Сеттона. Численный

алгоритм и его программная реализация тестировались путем сопоставления ре-

зультатов расчетов с результатами, отвечающими аналитическому решению.

Задача Сеттона [182] представляет собой задачу расчета характеристик ста-

ционарного процесса испарения жидкости при определенных допущениях на со-

стояние атмосферы. Граница жидкости ограничена координатой z = 0. Расчетная

область — полупространство x ≥ 0. Граница жидкости обдувается потоком воз-

духа, направленным вдоль оси x. Скорость воздушного потока — u(z). Положим,

что в атмосфере установилось некоторое стационарное распределение пара. При

этом у поверхности жидкости концентрация вещества равна некоторой равновес-

ной концентрации c

. Распределение концентрации пара в этом случае описывается

стационарным уравнением диффузии:

u(z)

∂c

∂x

∂

∂z

(z)

∂c

∂z

, (5.4.16)

в котором исключены горизонтальные турбулентные потоки в силу малости гори-

зонтальных градиентов концентрации по сравнению с вертикальными. Граничные

условия для этого уравнения имеют следующий вид:

c =

(

, при z = 0, x ≥ 0,

, при z → ∞, x ≥ 0,

(5.4.17)

где c

— концентрация в той части атмосферы, на которую граница жидкости не

оказывает влияния. Данная задача имеет аналитическое решение, если K

и u

являются степенными функциями z: u = az

, K

= bz

Введем безразмерную концентрацию χ =

c−c

−c

и автомодельную переменную

ξ =

b(2+m−n)

2+m−n

. Тогда решение может быть найдено в классе автомодельных

функций:

χ =

∞

−

m+1

2+m−n

−y

dy/Γ



1 −

m + 1

2 + m − n



, (5.4.18)

где Γ(x) — гамма-функция.

Определим константы a, b, n согласно [182]:

a =

5, 5u

∗

, b =

∗

5, 5m

, n = 1 − m, (5.4.19)

где динамическая скорость u

∗

= u

/(5, 5 ·2

) определена через параметр шеро-

ховатости и скорость ветра u

на высоте 2 м.

На рис. 5.10 представлены графики, отвечающие аналитическому решению за-

дачи Сеттона и решению с использованием представленной параллельной числен-

ной схемы по методу установления [176]. Результаты тестовых расчетов иллюст-

рируют хорошую сходимость численного решения к решению дифференциальной

задачи. В расчетах использована неравномерная по z разностная сетка с экспонен-

циальным сгущением в области больших градиентов концентрации, т. е. вблизи

подстилающей поверхности.

128 Компьютерная система прогноза последствий аварий

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0 500 1000 1500 2000 2500 3000 3500 4000

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0 500 1000 1500 2000 2500 3000 3500 4000

PSfrag replacements

численное решениечисленное решение

аналитическое решениеаналитическое решение

x, мx, м

c(x)/c

Рис. 5.10. Сопоставление аналитического и численного решения задачи Сеттона

= 2 м/с, z

= 0, 0213 см, m = 1/7, c

= 1, c

= 0) на высоте 2 м

5.5. Подсистема «Риск»

Подсистема предназначена для расчета аварийного риска по методологии, пред-

ставленной в главе 1, анализа его структуры и неопределенности. Реализация под-

системы допускает ее использование в виде самостоятельного продукта.

5.5.1. Архитектура. На рис. 5.11 проиллюстрирована архитектура подсисте-

мы «Риск», выделены два варианта работы системы, обусловленные различной

стратегией распараллеливания расчетов (см. раздел. 5.5.3).

В варианте 1 (архитектура «мастер»-«рабочий») на компьютере эксперта ра-

ботает интерпретатор языка Python в качестве управляющей программы, который

посредством динамически загружаемой библиотеки распределяет часть расчетов на

расчетные узлы. Вычислительно-затратные части анализа риска проводятся только

на расчетных узлах отдельной программой. Загружаемая библиотека и расчетная

программа реализованы на языке C++, что обеспечивает достаточно высокую вы-

числительную эффективность. Обмен данными между базовой библиотекой и рас-

четной программой осуществляется при помощи интерфейса посылки сообщений

MPI. Все операции ввода-вывода осуществляются на узле эксперта.

В варианте 2 узел эксперта служит лишь для распеределения заданий по се-

ти. На расчетных узлах работают интерпретаторы Python, однако все расчеты

проводятся при помощи загружаемой в Python базовой библиотеки. Файловая ба-

за данных (с рассчитанными полями ущерба) распределяется на расчетные узлы

при помощи любой файловой системы. Такая архитектура позволяет эффективно

5.5. Подсистема «Риск» 129

Компьютер эксперта

Расчетные компьютеры

База данных

Базовая библиотека,

C++

Расчетная

программа, C++

MPI

Модуль Python

Python

Вариант 1

Компьютер эксперта Расчетные компьютеры

Базовая библиотека,

C++

Модуль Python

Python

Сетевая

файловая система

База данных База данных

Вариант 2

Рис. 5.11. Архитектура подсистемы «Риск»

распараллелить решение задачи системными средствами и возможностями интер-

претируемого языка, не прибегая к дополнительному низкоуровневому программи-

рованию при помощи MPI.

5.5.2. Объектная модель. Модель подсистемы состоит из объектов двух се-

мейств классов. Базовое семейство классов отражает прагматику оценок риска.

Вспомогательное семейство позволяет представить аварийные ситуации в виде де-

ревьев, упрощающих анализ риска. Такое разделение на семейства обусловлено их

программной реализацией. Базовые операции при расчетах риска являются вычис-

лительно сложными и, следовательно, должны реализовываться на компилируемом

языке. Вспомогательное семейство, напротив, легче реализовать на бестиповом ин-

терпретируемом языке, позволяющем легко работать с такими типами данных, как,

например, деревья, словари и списки.

Базовое семейство классов. Анализ принятой методологии оценок аварийного

риска (см. раздел 1.2) позволяет выделить три класса основных объектов, которые

достаточно хорошо отражают суть прогноза риска (см. табл. 5.2).

Во-первых, при расчете условных оценок риска (1.2.4) используются различ-

ные распределения вероятностей: плотность распределения вероятности скорости и

направления ветра, температуры, состояния устойчивости атмосферы, вероятности

130 Компьютерная система прогноза последствий аварий

выпадения или невыпадения осадков. В подсистеме «Риск» указанные распределе-

ния представлены в виде объектов класса «Распределение».

Во-вторых, в условную оценку риска входят двумерные поля ущерба. При-

чем, сами модели ущерба можно рассматривать как фильтр (функцию), задан-

ный на двумерном поле поражающих факторов (токсических нагрузок, давления,

теплового потока и др.). Для удобства расчетов непрерывные поля поражающих

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

прямоугольной или полярной сетки. Абстракцией такой сеточной функции являет-

ся класс «Сетка». Третий класс объектов — абстрактный класс «Модель ущерба»,

является абстрактным интерфейсом для программного воплощения различных кон-

кретных моделей ущерба (см. раздел. 1.3), реализованных в виде дочерних классов

класса «Модель ущерба». Методами класса «Сетка» являются двумерные транс-

ляции и вращения, умножения на скаляр и поле, прибавление скаляра и поля,

интегрирование, фильтрация при помощи выбранной функции ущерба и др.

Все характеристики риска можно получить при помощи операций с объектами

трех вышеуказанных базовых классов. Очевидно, что такие опрерации являются

вычислительно-затратными. В целях вычислительной эффективности базовое се-

мейство классов было реализовано в виде библиотеки на языке C++ с организацией

интерфейса с языком Python.

Вспомогательное семейство классов. Семейство предназначено для пред-

ставления аварийных сценариев в виде логической упорядоченной структуры с це-

лью автоматизации и упрощения оценок риска, автоматизации анализа его струк-

туры.

Основные классы, участвующие в построении объектной модели реализаций

аварийной ситуации, показаны на рис. 5.12.

Для иллюстрации структуры риска рассмотрим условную аварию «дивер-

сия» на гипотетическом объекте утилизации химического оружия, состоящего из

нескольких производственных корпусов и подъездной железнодорожной ветки, рас-

положенной на охраняемой территории объекта.

Диверсия может реализоваться с определенной интенсивностью. Абстракци-

ей инцидента верхнего уровня — диверсии, является класс «Авария», имеющий

такие атрибуты, как, например, интенсивность и описание. Последствия от ди-

версии будут зависеть от того, на каком технологическом участке она произошла,

поскольку на различных участках боеприпасы находятся в различных количествах

и условиях, а, следовательно, параметры источника опасности (и, естественно, по-

следствия) будут различаться. Таким образом, аварию можно представить в виде

нескольких локализованных аварий (класс «Локальная авария»), каждая из ко-

торых имеет свою вероятность возникновения (не образующую полную группу в

общем случае). Предположим, что диверсия произойдет на участке железнодорож-

ного полотна (линейный участок). В связи с этим, места реализации локальной

аварии можно аппроксимировать несколькими точками с определенными весами

(класс «Точка»).

Кроме того, локальная авария может «проявится» несколькими способами (бое-

припасы какого типа попадают в зону воздействия), абстракцией которых является