CC BY
20
УДК 004.032.26
Моделирование процессов розлива нефтепродуктов по водной поверхности с использованием суперкомпьютеров с графическими ускорителями
И.А. Малый1, О.В. Потемкина1, С.Г. Сидоров2, А.С. Мочалов2, В.В. Неткачев2, И.Ф. Ясинский2, Ф.Н. Ясинский2 1 ФГБУВПО «Ивановский институт ГПС МЧС России», г. Иваново, Российская Федерация 2 ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина»,
г. Иваново, Российская Федерация
Авторское резюме
Состояние вопроса: Аварии с розливом нефти по водной поверхности, сопровождающиеся тяжелыми экологическими последствиями, происходят с досадной частотой. Существующие методы математического моделирования распространения нефтяного пятна требуют больших затрат машинного времени и отличаются невысокой точностью. Частично это можно объяснить высокой сложностью компьютерного моделирования нелинейных турбулентных процессов в трех измерениях.
Материалы и методы: Для решения соответствующих гидродинамических уравнений применяются численные методы расщепления, распараллеливания, неявные и полуявные схемы, скалярные прогонки. Для ускорения вычислений используется технология программирования графических процессоров CUDA. При реализации указанные численные методы модифицируются в соответствии с требованиями программной системы CUDA.
Результаты: Предложены алгоритмы и программы, позволяющие с существенным опережением и более высокой точностью предсказать движение и расширение нефтяного пятна, чем существующие программы, и высвободить время для проведения мероприятий, минимизирующих экологические последствия таких катастроф.
Выводы: Разработанная программа является частью системы, которая в настоящее время используется в МЧС в качестве тренажера и в перспективе будет внедряться в практику.
Ключевые слова: математическое моделирование, розлив нефти, параллельные вычисления, прогонка, расщепление, графические ускорители, CUDa.
Modelling of oil spill on a water surface by using supercomputers
with graphic accelerators
I.A. Maliy1, O.V. Potemkina1, S.G. Sidorov2, A.S. Mochalov2, V.V.Netkachev2, I.F.Yasinskiy2, F.N.Yasinskiy2 1 Ivanovo Institute of State Fire Service of Russian Ministry of Emergency Situations, Ivanovo, Russian Federation 2 Ivanovo State Power Engineering University, Ivanovo, Russian Federation
Abstract
Background: Accidents resulting in oil spill on water surface accompanied by serious environmental consequences are rather frequent. The current methods of mathematical modeling of oil spill spreading are computer time consuming and not precise enough. This can be attributed to the high complexity of 3-dimensional computer modeling of nonlinear turbulent processes.
Materials and methods: The corresponding hydrodynamic equations are solved by applying numerical techniques of splitting, paralleling, implicit and semi-implicit schemes, scalar marching. The computation speed is increased by the GPU programming technology CUDA. The above-mentioned numerical techniques are modified according to the CUDA programming requirements.
Results: The proposed algorithms and programs allow us to predict oil spill spreading much faster and more precisely than the current programs and to start taking measures reducing the environmental consequences of such accidents more quickly.
Conclusions: The developed program is part of the system which is currently used by the Ministry of Emergency Situations as a training simulator and will be put into practice in future.
Key words: mathematical modeling, oil spill, parallel calculations, marching, splitting, graphic accelerators, CUDA.
Существующие методы моделирования розлива нефти по водной поверхности не позволяют своевременно получить результат, необходимый для принятия оперативных решений по ликвидации последствий такой чрезвычайной ситуации. Поэтому необходима модель, которая будет давать адекватный результат, требовать минимальных затрат
машинного времени и хорошо распараллеливаться.
Нами были разработаны методы расчета и компьютерные программы, позволяющие: 1. При заданных конфигурации дна водоема, движении воды и ветра вычислить поле скоростей воды на поверхности водоема.
2. При заданном ветре и полученном поле скоростей воды на поверхности промоделировать движение нефтяного пятна, оценить его площадь и периметр, направление и скорость движения нефти, время достижения пятном определенных мест (например, берегов).
3. Определить потерю массы пролитых нефтепродуктов и изменение их состава за счет испарения более легких фракций.
4. Оценить уровень турбулентности в водоеме с помощью алгебраической и двух-параметрической дифференциальной моделей для расчета диффузии нефти из пятна в окружающую водную среду.
Исходные гидродинамические уравнения и методы их решения.
1. Вычисление полей скорости в воде выполнено с помощью расширенного метода полного потока Штокмана [1]. Расширение состоит в добавлении нами в динамические уравнения конвективных членов:
т13
ди ди ди , дк —1 + и1 —1 + и2 —1 - 1и2 = д —-
дх1
' дх2
дх1 рН
(я 2,
V
д2и д2ы.
2, ^
(1)
ди2
~дГ
ди2
'"дхТ
дх{ дх2 _
ди2 . дк - и2 —2 - 1и1 = д —-2 дх2 дх2
т23
рН
-Ни,
дк д:
= Н
дх2 ди
д 2и2 дх|
дх1 дх2
дН_ дхЛ
дН
дх,
(2)
(3)
(4)
Т13 = р( -и1); Т23 = Р(( -и2). (5)
Здесь : - время; х1, х2 - декартовы координаты на поверхности воды; и1, и2 - составляющие осредненной по глубине скорости течения; №2 - составляющие скорости ветра; Н(х1,х2) - расстояние от поверхности до дна водоема; т13, т23 - составляющие ветрового напряжения на поверхности воды; к(:,х1,х2) -уровень водной поверхности; А, уЭ, V - кинематические вязкости; д - ускорение свободного падения; I - константа Кориолиса; Н, р -постоянные.
Решение ищется с помощью расщепления по процессам и направлениям (в пределах малого шага по времени). Например, вместо (1) получим следующую цепочку уравнений:
ди1 ~дТ ди1 + и1—1 1 дх1 = А д 2и21; дх( (6)
ди1 ~дТ ди1 + и2—1 2 дх2 = А д 2и21; дх| (7)
ди1 ~дТ дк =д^- дх1 + !13-Хэ Ни1 - 1и2. рН V 12 (8)
При этом конвективные слагаемые заменяются разностными приближениями с использованием «противоточных производных»:
ди1
'"дх
. " Ш ((
1 у" + |и1/у
,к+1 >1ц
- и.
к+1 ^ 11-1]
,,к+1 и и1+1 ] - и1]
] )).
(9)
+ ( -\и1ц
Здесь , ] - горизонтальная и вертикальная нумерация узлов сетки, наброшенной на рассматриваемую область; верхний индекс - номер момента времени (расстановкой этого индекса линеаризуются конвективные слагаемые); Л - шаг квадратной сетки.
Составляющие градиента аппроксимируются следующим образом: дк
дк = 2Л (к+1]+1 + к+1]-1 - к'-1]+1 - к-1]-1)'
(10)
Разностные выражения приводятся к
виду
к,,к , пк,,к+1 , к,,к+1 як /лл\
а1 ¡]и1-1 у + Р и1 ¡] + У1 ]и1+1 у = 8 ¡] , (11)
к которому применяется скалярная прогонка.
На прямом ходе вычисляются прогоноч-ные коэффициенты:
4+1У =-Ук / (вкк + 44); (12)
Чк+1] = (8к -аумк )/( + ). (13)
На обратном ходе находятся искомые величины:
:. (14)
Аналогичные операции выполняются в направлении х2.
2. Движение нефтяного пятна подчиняется уравнениям:
+ц дЧ + ц2 V = д: 1 2
цк+1 = ¡к/¡к и1-1 ] = ¡-¡]и1]
дх1 дх2
= -д
■(к + А)-р1 (ц - и1) + й ^
+ц +ц2 =
д:
= -д
дхл
дх,
дх-
-(к + А)-р1 ( - и2 ) + й
дг2
д У2
г дг2
д:
= А
-А + V
(я 2,
дх2
д ц2
-А + д А
дх1 дх2
+ д2А^ й ^ г дг2
"+дх| , +
(15)
(16)
(17)
- Ои + О+
Здесь Ц1, V - скорости нефтяного пятна; А -толщина нефтяного пятна; Ои - скорость изменения толщины пятна за счет его испарения; О+ - приток нефти за счет подводного источника, если таковой имеется; слагаемые р1(^1 - и1), р1(Ц2 - и2) - учитывают трение нефтяного пятна при его относительном движении по поверхности воды.
Решение этой системы выполнялось аналогично изложенному выше.
Для турбулентной кинематической вязкости справедлива оценка
й2 = 0,8 Ь2/ъ , (18)
где Ь - энергия турбулентных пульсаций; е -скорость диссипации турбулентности. Они подчиняются уравнениям:
- = Dz dt z
dV1 dz
^ = 1,38 • Dz
dt z
dV1 dz
dV2 I21 д
dz J T Hz
l2 + ^5V2 J21
J + v dz J
db_ ~dz
(19)
(20)
d
D
db
-1,4-
д7 { 7 д7) Ь и граничным условиям (ось I направлена с поверхности воды вертикально вниз): • при I = 0
0 5Ц = _Тз; 0 У = т23 ;
2 д7 р0 ' 2 81 '
Ро
(23)
V2 = 0,03 • W2 +1,1 • u2.
Наши численные эксперименты подтверждают это.
Также достаточно хорошей оценкой для кинематической турбулентной вязкости, является величина, полученная с помощью метода Л. Прандтля, «О длине пути перемешивания»:
Dz = 0,16^V2 + V22 . (24)
3. Учет потери части массы нефтяного пятна за счет испарения легких фракций выполнен на основании работы2 и материалов, предоставленных фирмой ОАО «Транзас Морские Технологии»3.
Если известны градусы API нефти [2], то весь объем пролитой нефти делится на пять равных по объему частей. Для них вычисляются:
Материалы предоставлены фирмой ОАО «Транзас Морские Технологии» при личном общении.
2 NOAA, 2000. ADIOS™ (Automated Data Inquiry for Oil Spills) version 2.0. Seattle: Hazardous Materials Response and Assessment Division, NJAA. Prepared for the U.S. Coast Guard Research and Development Center, Groton Connecticut.
3 Материалы предоставлены фирмой ОАО «Транзас Мор-
ские Технологии» при личном общении.
• температуры кипения
Г = 457,16 - 3,3447 • API + (0,2 • i - 0,1)
х (1356,7 - 247,36 • ln( API)), i = 1,2,3,4,5;
• молярный объем
Ni = 7 • 10-5 - 2,102 • 10-7 • Г1 +10-9 • Гf;
• давление пара
P = 103 • expI-(4,4 +1,45• lnГ,)х
1,8031 Гт- 1j-0,803 • lniГ-
(26)
(27)
где Т - температура воды.
Испарение существенно зависит от ско рости ветра:
[0,005 при М < 4 ,
C(W) = <
0,002 • W°,78 при W > 4 м
(28)
й7 — = 0; й7 — = 0; (21)
7 д7 7 д7 -
• при I = Н
V = 0; У2 = 0; Ь = 0; е = 0. (22)
Турбулентные уравнения решались до установления с помощью явной схемы с шагом по времени т1, резко меньшим, чем шаг т для систем (1)-(5), (15)-(17).
Отметим, что, согласно опыту фирмы «Транзас Морские Технологии»1, оказалась весьма удобной следующая зависимость, найденная эмпирически: У = 0,03 • М1, +1,1 • и1,
Скорость уменьшения массы каждой из п = 5 частей описывается уравнением
РСРт / = 1,2.....5, (29)
* РТУ5
=1 р ¡ы,
где т, - масса каждой из пяти частей; Р -площадь нефтяного пятна; р, - плотность / части; Р - газовая постоянная.
Интегрирование этой системы выполняется элементарно.
Отметим, что, согласно [2], под действием ветра нефтяное пятно от единичного кратковременного источника принимает форму эллипса:
• с полуосями
A = 1,417
Pw -Р0 Pw
Уб f V У3 tyA.
V 3/
A2 = A, + 0,15• W/3 t/4;
• площадью S = n • A,A2;
• периметром
p = n {o,15 (A, + A2) -yjA,A2}.
Здесь vw - кинематическая вязкость воды; V -объем пролитой нефти; t - время.
Программная реализация модели розлива нефти. Изложенные вычислительные схемы были реализованы в виде программы на языке C++ с использованием технологии NVIDIA CUDA.
Для вычисления противоточных производных был использован прямой подход доступа к разделяемой памяти [9]. Эксперименты показывают, что такой подход обеспечивает максимальную скорость расчета явных схем, по сравнению с другими подходами.
Для вычисления скоростей воды и нефти с помощью прогонки применен подход раз-
деления по линиям [10]. Данный подход в чистом виде (при вычислении только одной величины) дает прирост производительности примерно в 4 раза (рис. 1). Еще одним преимуществом данного подхода является небольшое количество блоков сетки и потоков блоков программной модели CUDA, что позволяет запускать расчет больших сеток даже на слабых графических процессорах, в которых максимальное количество потоков ограничено.
7 -
Рис. 1. Производительность прогонки с разделением по линиям в чистом виде (при вычислении одной величины)
Над программой были произведены нагрузочные тесты на видеокарте NVIDIA GeForce 8600M, по результатам которых можно сделать вывод, что скорость вычисления, по сравнению с последовательным вариантом, увеличилась примерно в 5 раз (рис. 2).
Стоит отметить, что видеокарта NVIDIA GeForce 8600M является одной из самых слабых из тех, которые поддерживают технологию CUDA. На более мощных видеокартах ожидается еще большее ускорение. Соответствующие эксперименты планируется провести в перспективе.
^<ГЧСП^Г1ЛСОГ-СОСПС=1С=1С=1С=1С=1С=1С=1С=1С=1С=1С=1
Рис. 2. Ускорение вычислений за счет применения технологии CUDA
Для реализованной программы была разработана графическая оболочка в виде ГИС. В качестве ГИС движка используется DotSpatial. Данная система разработана на платформе .NET. Расчетный модуль интегрируется с ней с помощью технологии CLI. Эксперименты, проведенные в этой программе, показывают адекватность математической модели: скорость, с которой происходит распространение нефти, соответствует результатам реальных наблюдений, полученным у компании «Транзас Морские Технологии»4, у которой есть опыт в наблюдении и моделировании розлива нефти по водной поверхности.
В данной ГИС есть возможность графического ввода исходного нефтяного пятна в виде изображения и последующего вывода в каждый следующих момент времени также в виде изображения. Это особенно удобно, если слежение за развитием чрезвычайной ситуации ведется с помощью беспилотных летательных аппаратов.
На рис. 3 представлены последовательные стадии расширения нефтяных пятен (обозначены черным цветом) по водной поверхности около трех островов (обозначены светлым). Два нефтяных пятна сливаются в одно.
\ '1 1 J . \ X X "X
\ ■ < 1 i . \ ч -х X "К
! i , \ X "V X "К
Рис. 3. Моделирование розлива нефти по водной поверхности в ГИС
Также в этой системе есть возможность выгрузки данных в формат KML, который поддерживает компания Google. Это дает возможность наблюдать за распространением
4 Материалы предоставлены фирмой ОАО «Транзас Морские Технологии» при личном общении
нефти в других ГИС. Например, динамика распространения нефти в программе Google Earth представлена на рис. 3.
Рис. 3. Динамика розлива нефти по водной поверхности в программе Google Earth
Заключение
Использование графических ускорителей системы CUDA, а также предложенных алгоритмов позволяет с минимальной затратой времени предсказать последствия розлива нефти по водной поверхности (в том числе, время достижения ею берегов водоема или других критически важных объектов).
Таким образом, сотрудники МЧС получают необходимое время для принятия решений по ликвидации последствий данной чрезвычайной ситуации.
Список литературы
1. Кочергин В.П. Теория и методы расчета океанических течений. - М.: Наука, 1978. - 126 с.
2. Fay J.A. Physical processes in the spread of oil on a water surface. Proc. On Prevention and Control of Oil Spill / American Petroleum Institute. - Washington, DC, 1971. -Р. 463-467.
3. Brebbia C.A. Oil Spill Modelling and Processes // WIT Press. - United Kingdom, 2001.
4. Tkalin A.V. Evaporation of petroleum hydrocarbons from films on a smooth seasurface // Oceanology of the Academy of Sciences of the USSR. - 1986. - Vol. 26. -Р. 473-474.
5. Stiver W., Mackay D. Evaporation rate of spills hydrocarbons and petroleum mixtures // Environmental Science and Technology. - 1984. - Vol. 18. - Р. 834-840.
6. Mackay D., Buistt I.A., Marcarenhas R., Pater-son S. Oil spill processes and models: Environment Canada Manuscript Report. - Ottawa: Ontario, 1980. - No. EE-8.
7. Mooney M. The viscosity of a concentrated suspension of spherical particles // J. Collodal Science. - 1951. -Vol. 10. - Р. 162-170.
8. Garo J.P., Vantelon J.P., Gandhi S., Torero J.L. Some observation on the preboilover burning of a slick of oil on water // Proceeding of the Nineteenth Arctic and Marine Oil Spill Program (AMOP) Technical Seminar. - Ottawa, 1996. - Р. 1611-1626.
9. Неткачев В.В. Технология CUDA и ее применение для решения задач: учеб. пособие. - Иваново, 2013.
10. Неткачев В.В. Моделирование механики жидкости и газа с применением технологии CUDA // Вестник ИГЭУ. - 2013. - Вып. 3. - С. 69-73.
References
1. Kochergin, V.P. Teoriya i metody rascheta okeani-cheskikh techeniy [Theory and methods of calculating ocean currents]. Moscow, Nauka, 1978. 126 p.
2. Fay, J.A. Physical processes in the spread of oil on a water surface. Proc. On Prevention and Control of Oil Spill. American Petroleum Institute. Washington, DC, 1971, pp. 463-467.
3. Brebbia, C.A. Oil Spill Modelling and Processes. WIT Press. United Kingdom, 2001.
4. Tkalin, A.V. Evaporation of petroleum hydrocarbons from films on a smooth seasurface. Oceanology of the Academy of Sciences of the USSR, 1986, vol. 26, pp. 473-474.
5. Stiver, W., Mackay, D. Evaporation rate of spills hydrocarbons and petroleum mixtures. Environmental Science and Technology, 1984, vol. 18, pp. 834-840.
6. Mackay, D., Buistt, I.A., Marcarenhas, R., Pater-son, S. Oil spill processes and models: Environment Canada Manuscript Report. Ottawa, Ontario, 1980, no. EE-8.
7. Mooney, M. The viscosity of a concentrated suspension of spherical particles. J. Collodal Science, 1951, vol. 10, pp. 162-170.
8. Garo, J.P., Vantelon, J.P., Gandhi, S., Torero, J.L. Some observation on the preboilover burning of a slick of oil on water. Proceeding of the Nineteenth Arctic and Marine Oil Spill Program (AMOP) Technical Seminar. Ottawa, 1996, pp. 1611-1626.
9. Netkachev, V.V. Tekhnologiya CUDA i ee prime-nenie dlya resheniya zadach [CUDA technology and its application to solving problems]. Ivanovo, 2013.
10. Netkachev, V.V. Modelirovanie mekhaniki zhid-kosti i gaza s primeneniem tekhnologii CUDA [Modelling of fluid mechanics by CUDA technology]. Vestnik IGEU, 2013, issue 3, pp. 69-73.
Малый Игорь Александрович,
Ивановский институт государственной противопожарной службы МЧС России, начальник института, генерал-майор внутренней службы, e-mail: [email protected]
Потёмкина Ольга Владимировна,
Ивановский институт государственной противопожарной службы МЧС России, заместитель начальника института по учебной работе, майор внутренней службы, e-mail: [email protected]
Сидоров Сергей Георгиевич,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», кандидат технических наук, доцент, зав. кафедрой высокопроизводительных вычислительных систем, e-mail: [email protected]
Мочалов Александр Сергеевич,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», программист отдела геоинформационных технологий ИВЦ, e-mail: [email protected]
Неткачев Владимир Владимирович,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», преподаватель кафедры высокопроизводительных вычислительных систем, ведущий программист, e-mail: [email protected]
Ясинский Игорь Федорович,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», преподаватель кафедры высокопроизводительных вычислительных систем, e-mail: [email protected]
Ясинский Федор Николаевич,
ФГБОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина», профессор кафедры высокопроизводительных вычислительных систем, телефон (4932) 29-98-29.
