Имитационная модель распространения поражающего агента в городской застройке Текст научной статьи по специальности «Физика»

ИМИТАЦИОННАЯ МОДЕЛЬ РАСПРОСТРАНЕНИЯ ПОРАЖАЮЩЕГО АГЕНТА В ГОРОДСКОЙ ЗАСТРОЙКЕ

Самохина А. С.

(Институт проблем управления им. В.А. Трапезникова РАН, Москва) absamokhin@yandex.ru Сетуха А. В., Кирякин В. Ю. (Военно-воздушная инженерная академия им. проф. Жуковского, Москва) setuhaav@rumbler.ru, Kiryakin@rambler.ru Марчевский И. К., Щеглов Г. А.

(Московский государственный технический университет им. Н. Э. Баумана, Москва) iliamarchevsky@mail.ru, georg@energomen.ru

Разработаны теоретические основы и создано программное обеспечение компьютерного объемного моделирования распространения биологических поражающих агентов импульсного выброса в условиях сложной городской застройки с учетом неравномерности ветрового потока, возникающего при обтекании зданий. Рассмотрено применение экспертной системы для определения и согласования информативных признаков входных данных численного моделирования и классификации опасных зон распространения.

Ключевые слова: экспертные системы, компьютерное объемное моделирование, биологический агент.

Введение

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

телей инфекционных болезней и токсинов. Ничтожно малые инфицирующие дозы, отсутствие высокочувствительных и специфичных методов и средств экспресс-индикации микроорганизмов в пробах окружающей среды, недостаточная эффективность препаратов общей экстренной профилактики и этиопато-генетического лечения определяют потенциальную угрозу использования биологических поражающих агентов (БПА) для нарушения безопасности страны и создания чрезвычайной ситуации, вызванной биологическими факторами (БЧС) [6, 7]. В этой связи создание системы принятия решений (СПР) для предупреждения, противодействия и ликвидации последствий БЧС и разработка ее информационного и компьютерного обеспечения является актуальной и своевременной проблемой [2].

Прогнозирование санитарно-эпидемиологической обстановки является обязательным элементом создаваемой СПР, включающей систему подготовки данных для планирования оптимальных профилактических и противоэпидемических мероприятий в целях минимизации возможного ущерба и ликвидации медико-санитарных последствий в кратчайший период. Одной из важнейших задач прогнозирования является определение границы зоны биологического заражения и связанная с ней задача определения уровня инфекционной заболеваемости среди населения, попавшего в эту зону, определения зараженных поверхностей и площадей, а также подсчет сил и средств, необходимых для проведения дезинфекции зараженной территории [2].

При построении имитационной модели распространения аэрозольного или мелко-дисперсионного пылевого облака при условии импульсного выброса мелко-дисперсионного биологического агента заражения в многоэтажной жилой застройке эксперты играют ведущую роль при выборе и согласовании информативных признаков входных параметров модели, определении степени зараженности территорий, классификации опасных зон и выборе способов дезинфекции после определения системой поверхностей и площадей, подлежащих дезинфекции.

Входными параметрами модели являются: место проведения выброса, которое определяется на основании сведений, полученных санитарно-эпидемической разведкой; диаметр и масса частиц БПА, его структура (аэрозоль или пыль); метеорологиче-

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

1. Общая постановка задачи распространения мелкодисперсионного биологического агента

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

Предполагается, что группа зданий обтекается ветровым потоком, имеющим заданную постоянную скорость и направление на большом удалении от зданий. Локальное поле скоростей в непосредственной близости от зданий определяется на основе прямого численного моделирования. Считается, что движение частиц не возмущает это поле скоростей. Предполагается, что частицы примеси «рождаются» в заданной области пространства и движутся под воздействием силы тяжести и аэродинамических сил, обеспечивающих их перенос в воздушном потоке. Расчет движения частиц производится до момента их соприкосновения с землей или поверхностью зданий.

Для проведения расчетов и анализа их результатов вводится неподвижная система координат, в которой ось Oz направлена по нормали к поверхности земли. Траектории движения частиц примеси определяются внутри расчетной области, которая задается неравенствами хх < х < х2, у < у < у2, 0 < z < Н. Нижняя грань расчетной области z = 0 считается непроницаемым экраном, который моделирует поверхность земли.

Рис. 1. Расчетная область

В результате определяется концентрация осажденной примеси на земле и на поверхностях зданий и сооружений.

Кроме того, в области, задающей место «рождения» частиц, выделяются зоны, указывающие местоположение тех «родившихся» частиц, которые в процессе дальнейшего движения осаждаются на землю, здания или уносятся из расчетной области.

При расчете области заражения вначале осуществляется ввод геометрии объектов, расположенных на рассматриваемом участке местности, после этого вводится скорость и направление ветра, представленные экспертами, а также согласованные параметры БПА. Результатом данного этапа работы является построение трехмерной модели зданий с распределением вектора скоростей аэрозольного или пылевого облака агента и границами аэрационных зон [4].

Допущение об отсутствии влияния примесей на поле скоростей воздушного потока позволяет разделить процесс решения задачи на два этапа:

1) расчет поля скоростей ветрового потока, возникающего при обтекании зданий и сооружений;

2) расчет движения примеси в уже рассчитанном поле скоростей.

Отметим, что возникающее поле скоростей является, вообще говоря, нестационарным. В связи с этим возможно несколько вариантов решения общей задачи, при которых для нахождения траекторий движения частиц используется:

а) нестационарное поле скоростей, полученное на определенном временном интервале, соответствующем временному интервалу движения частиц;

б) стационарное поле скоростей, соответствующее определенному моменту времени при решении аэродинамической задачи;

в) стационарное поле скоростей, получаемое осреднением нестационарного поля скоростей по времени;

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

Заметим, что вариант а) требует значительных вычислительных затрат при моделировании движения частиц на втором этапе. Варианты б) и в) позволяют применить методы ускоренного расчета траекторий движения частиц, которые будут описаны ниже, и за счет этого существенно сократить вычислительные затраты. Вариант г) является компромиссным вариантом, требующим существенно меньше вычислительного времени для расчета траекторий частиц, чем при варианте а), и при этом позволяет внести элемент нестационарности в поле скоростей ветрового потока.

2. Метод решения аэродинамической задачи

Поле скоростей ветрового потока находится путем прямого численного моделирования процесса обтекания группы зданий на основе вихревой модели идеальной несжимаемой жидкости. При этом использовался метод вихревых рамок, описанный в [4]. Предполагается, что на здания набегает воздушный поток, имеющий постоянную скорость на бесконечности. На поверхностях тел ставится условие непротекания (условие равенства ну-

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

В используемом методе поверхности обтекаемых тел и вихревые следы заменяются дискретными вихревыми элементами. Интенсивности вихревых элементов (вихревых рамок), аппроксимирующих на каждом шаге интегрирования по времени поверхности обтекаемых тел, находятся из условия непротекания. Движение вихревых элементов, сошедших в поток с линий отрыва, описывается на основе закона о движении вихревых линий с жидкими частицами и о сохранении циркуляции скорости по контуру, движущемуся вместе с жидкостью.

В результате расчетов на каждом шаге интегрирования по времени определяются интенсивности вихревых рамок, моделирующих поверхности обтекаемых тел, а также интенсивности и законы движения вихревых элементов, моделирующих вихревой след. С использованием закона Био-Савара можно определить значение вектора скорости в любой точке течения жидкости.

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

3. Определение траектории движения примеси биологического агента в рассчитанном поле скоростей

Начальное положение частиц примеси в расчетной области предполагается заданным, а начальная скорость совпадает со скоростью потока в соответствующей точке. Движение частицы примеси в потоке среды с плотностью р0 и полем скоростей W ( г ) описывается функцией г ^), задающей закон движения частицы, и являющейся решением задачи Коши [1, 5]:

(1) тг = -2-CDdpo ^ -т7) ...г(0) = ге ,...г(0) = W(ге),

где т - масса частицы, й - диаметр частицы, Сг - коэффициент сопротивления частицы, ге - начальное положение частицы.

24

Считается, что Сг = — (формула Стокса), а число Рейнольдса Яе

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

й ■ # -г| . ..

диаметру, Яе = —--------, г и г - первая и вторая производные

V

функции г (/) по времени (скорость и ускорение частицы).

Осаждение примеси происходит в момент соприкосновения с поверхностью земли или здания.

В результате первого (аэродинамического) этапа решения задачи для каждого шага по времени находятся циркуляции вихревых рамок, моделирующих поверхности зданий и вихревой след, и положения вихревых рамок, моделирующих вихревой след. Поле скоростей определяется формулой - - N _

# (Г) = Гх+Т V (Г)

I=1

где г - радиус вектор точки пространства, в которой вычисляется вектор скорости воздуха # , - скорость невозмущенно-

го потока на бесконечности, N - общее количество вихревых рамок, моделирующих поверхности тел и вихревой след, V (г) -скорость, индуцируемая в рассматриваемой точке вихревой рамкой с номером г в соответствии с законом Био-Савара.

Для однократного вычисления скорости требуется вычисление влияния от N вихревых рамок, общее количество операций умножения и деления при этом порядка N. Для интегрирования движения Р частиц в течение S временных шагов потребуется выполнить порядка №Р-8 операций. Количество операций можно существенно уменьшить, если предварительно построить в расчетной области сетку с числом узлов О, в узлах которой предварительно вычислить скорости потока, а затем при помощи линейной интерполяции вычислять скорости в любой точке. При этом на вычисление скоростей в узлах потребуется порядка N0 операций, а на вычисление правых частей уравнения путем интерполяции - порядка Р8 операций.

Итак, при непосредственном вычислении скоростей требуется порядка №Р-8 операций, а с использованием сетки порядка

(N•0 + Р-5) операций. Видно, что при достаточно больших значениях параметра N использование сетки оправдано при условии, что О ^ Р-5.

В реальном расчете было использовано О ~ 104, Р ~ 105,

5 ~ 104, N ~ 104 N ~ 104 . Как видно, вычислительные затраты, при использовании сетки являются несоизмеримо малыми, по сравнению с прямыми вычислениями.

Интегрирование уравнений (1) производится методом Эйлера с коррекцией 2-го порядка точности. При этом имеются 4 условия остановки счета:

1) осаждение частицы на экран г = 0;

2) осаждение частицы на поверхность здания;

3) выход частицы за границу расчетной области;

4) превышение некоторого заданного количества шагов по времени.

Задача о расчете движения примеси допускает эффективное распараллеливание, поскольку вычисление скорости для каждого узла сетки может выполняться независимо, и интегрирование уравнения движения отдельной частицы также может выполняться независимо от остальных. Эффект от распараллеливания вычислительной процедуры можно оценить по формуле

1 т ■ 1п(п) О + Р

К к- + -

п а N0 + БР ’ где п - количество однотипных компьютеров в сети, а - характерное время выполнения одной арифметической операции; т -характерное время пересылки данных между двумя компьютерами.

Множитель 1п п определяется архитектурой используемой стандартной библиотеки МР1.

тз О + Р 4

В реальных расчетах множитель-------------10 .

N0 + БР

Алгоритм параллельного решения задачи включает в себя следующие стадии:

1) ввод исходных данных и их рассылка на все компьютеры

сети;

2) вычисление скоростей в узлах сетки на сети компьютеров (каждый компьютер вычисляет скорости в определенном подмножестве узлов сетки);

3) обмен результатами и рассылка на все компьютеры значений скорости во всех узлах сетки;

4) интегрирование уравнений движения частиц (каждый компьютер вычисляет траектории определенного подмножества частиц);

5) сбор и обработка результатов на «головном» компьютера

Результаты расчета представляют собой координаты частиц на момент окончания счета, а также информацию о том, какое из вышеуказанных условий остановки расчета реализовалось. По этим результатам вычисляется:

• Концентрация «на земле», которая определяется путем построения на экране прямоугольной сетки и определения числа частиц, попавших в каждую ячейку этой сетки. В качестве концентрации примеси будем вычислять отношение количества частиц в ячейке к общему количеству частиц, осевших на экране.

• Концентрация «на здании», которая определяется аналогичным образом, при этом все поверхности зданий разбиваются на прямоугольные ячейки.

• Диаграмма осаждения. Для качественного анализа распространения примеси представляет интерес зависимость результата движения частицы (какое из условий останова реализовалось) от начального положения частицы ге . Например, если точки старта частиц лежат в некоторой «контрольной» плоскости а, на этой плоскости могут быть выделены зоны, частицы из которых оседают на экране, здании или уходят из расчетной области.

4. Компьютерное моделирование распространения биологического агента в городской застройке

Расчётный пример. В области, показанной на рис. 2, размещены три здания: 1 - высотный дом, 2 - угловой дом, 3 - объ-

ект социальной инфраструктуры, например, детский сад. Отметим, что такое сочетание и взаиморасположение зданий является достаточно распространенным элементом городской застройки. В расчете используется источник частиц - эмиттер, представляющий собой прямоугольный параллелепипед с центром в точке х = -130; у = 50; г = 10. Длины сторон параллелепипеда равны соответственно Ах = 32; Ду = 100; Аг = 18. Эмиттер равномерно заполнен 40300 точками. Параметры частиц т = 2-10-12;

й = 10-4. Скорость набегающего потока равна Уж = 1. Распараллеливание вычислительной процедуры проводилось на 10 компьютерах.

На основании проведенного модельного расчета можно сделать следующие основные выводы:

1. Кроме фасада высотного дома, стоящего поперек потока, сильно загрязненными оказались углы зданий.

2. Неочевидным фактом является загрязнение теневой стены здания 3 (детский сад), поскольку на нее попадают частицы, закрученные вихревым потоком, образующимся за угловым зданием. Это показывает важность учета вихревых структур при расчете переноса загрязнений.

o'

Ш

ЭМИТТЕР

250 150 0 -50 -150

Рис. 2. Моделируемая застройка, расчетная область и эмиттер частиц

Рис. 3. Поле скоростей ветрового потока на уровне земли

К <0.4 0.4<7<1.1 1.1£Й<1.3 Р ^ 1.3 У=ПУШ

Рис. 4. Аэрационные зоны на уровне земли (см. табл. 1)

Г

Рис. 5. Диаграмма распределения модуля скорости ветрового потока на уровне земли

Рис. 6. Характерные траектории частиц. Изометрия с подветренной стороны

Рис. 7. Характерные траектории частиц. Вид сверху

Рис. 8. Характерные траектории частиц. Вид с наветренной

стороны

Рис. 9. Концентрации частиц на экране

КОНЦЕНТРАЦИЯ

Рис. 10. Концентрации на фасадах зданий. Вид с подветренной

стороны

КОНЦЕНТРАЦИЯ И ВЫСОКАЯ

Рис. 11. Концентрации на фасадах зданий. Вид с наветренной стороны

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

2. осаждение на стенах здания

3. выход из расчетной области без загрязнения

4. осаждение на поверхности экрана 7

150 0 -50

Рис. 12. Диаграмма осаждения. Ближний к зданиям слой эмиттера (Х = -114)

5. Применение системы принятия решений для определения поверхностей, подлежащих дезинфекции

СПР на основании полученных вычислительным модулем данных о загрязнении поверхностей зданий, сооружений и почвы, позволяет выделять разные классы аэрационных зон: застоя, комфортных условий, умеренного ускорения воздушного потока, дискомфортные зоны, что помогает принимать решение о дезинфекции территорий. Алгоритм работы СПР заключается в следующем.

На первом шаге эксперты при помощи компьютерной процедуры голосования согласовывают диапазоны значений коэффициента модуля скорости V = V / , где V - величина

скорости ветрового потока в рассматриваемой точке, VIX - скорость невозмущенного потока на бесконечности, которые используются для выделения аэрационных зон (зон значений модуля скорости). В нашем случае используется таблица 1 [4].

Таблица 1. Диапазоны значений коэффициента модуля скорости для аэрационных зон_____________________________________

Название аэрационной зо- Диапазон значений коэффициента

ны модуля скорости

Зона застоя 0 < V < 0.4

Нормальная зона 0.4 < V < 1.1

Допустимая зона 1.1 < V < 1.3

Дискомфортная зона V > 1.3

На втором шаге вычисляются поля скоростей ветрового потока в интересующих сечениях, которые отображаются в виде векторов, из рисунка которых получаем наглядную качественную картину воздушного потока в выбранном сечении (рис. 3-5). На мониторе компьютера высвечиваются цветовые диаграммы распределения модуля скорости в выбранном сечении (рис. 6-8), отображаются границы осаждения биологического агента (рис. 9-12).

На третьем шаге экспертами принимается решение о зонах дезинфекции. Для этого на экранах компьютера, на основании рис. 9-12, каждый эксперт выделяет зоны, которые он считает подлежащими дезинфекции. При этом он основывается на информации о характере агента, его живучести, индексе контаги-озности, метеорологических данных и т. д. Компьютерная система складывает все рисунки, полученные от экспертов, и выдает суммарный рисунок, на котором наиболее ярко окрашена зона безусловной дезинфекции, и далее по градации цвета зоны, набравшие соответствующее количество голосов экспертов. На основании информации о зонах, подлежащих дезинфекции, решается вопрос о силах и средствах, используемых при дезинфекции.

Заключение

Разработанные и реализованные методы вычисления распространения биологического поражающего агента в городской застройке с учетом вихревых потоков являются важными в системе предупреждения, преодоления и ликвидации последствий чрезвычайных ситуаций, обусловленных биологическими факторами. На реализованной компьютерной системе проводится моделирование развития чрезвычайной ситуации в городе при импульсном выбросе мелко-дисперсионного поражающего агента и принятие решений по ликвидации последствий чрезвычайных ситуаций с учетом мнений экспертов.

Литература

1. ВОЛКОВ К. Н. Разностные схемы интегрирования уравнений движения пробной частицы в потоке жидкости или газа // Вычислительные методы и программирование. М.: Изд-во Московского университета. 2004. Т. 5. № 1. С. 5-21.

2. ГАБУСУ П. А., МИХЕЕВ О. В., САМОХИНА А. С. Экспериментальные исследования имитационной математической модели распространения опасного эпидемического заболевания / Тр. междунар. конф. «Параллельные вычисле-

ния и задачи управления» РАСО 2006, Ин-т прбл. упр. - М., 2006. С. 1494-1508.

3. ГАБУСУ П. А., МИХЕЕВ О. В., САМОХИНА А. С. Идентификация сценариев эпидемического заболевания. // Проблемы безопасности и чрезвычайных ситуаций. 2006. № 2. С. 92106.

4. ГУТНИКОВ В. А., ЛИФАНОВ И. К., СЕТУХА А. В. О моделировании аэродинамики зданий и сооружений методом замкнутых вихревых рамок. // Механика жидкости и газа. 2006. №4. С. 79-93.

5. ЛОГАЧЕВ И. Н., ЛОГАЧЕВ К. И. Аэродинамические основы аспирации / СПб.: Химиздат. 2005. - 659С.

6. Организация и проведение противоэпидемических мероприятий при террористических актах с применением биологических агентов. МР 2510/11646-01-34.

7. Организация и проведение режимно-ограничительных мероприятий в зонах стихийных бедствий и техногенных катастроф. МУ 1.2.793-99.

Статья представлена к публикации членом редакционной коллегии В.В. Кульбой