CC BY
47
УДК 531.383:532.516 ИБО 65M06, 76Б05
Математическое и компьютерное моделирование эволюции
нефтяного пятна
А. В. Трепачева, Ф. Б. Буртыка
Южно-Российский 'региональный центр информатизации Южный федеральный университет пр. Стачки, д.200/1, корп.2, Ростов-на-Дону, Россия, 344090
В работе рассматривается модель переноса нефтяного пятна по поверхности воды, базирующаяся на квазилинейном уравнении конвекции-диффузии, составленном относительно толщины пятна. Для решения последнего используется метод конечных разностей. Двумерная задача сводится к одномерной с помощью метода переменных направлений. Для одномерного случая рассматриваются традиционные разностные аппроксимации, а также разностные схемы, которые можно получить с помощью метода автоматической генерации разностных схем, основанного на методах компьютерной алгебры. Для сгенерированных схем оценивается порядок аппроксимации и схемная вязкость. Для явных схем рассматривается вопрос об условиях устойчивости.
В случае отсутствия конвекции проводится численное сравнение автомодельного решения квазилинейного уравнения теплопроводности с традиционной неявной разностной схемы и одной неявной схемой, сгенерированной автоматически. Сравнение показало, что сгенерированная неявная схема позволяет получить относительную погрешность решения меньшую, чем для традиционной схемы.
В соответствии с полученной неявной схемой проводится расчёт изменения толщины пятна при условии наличия испарения нефти с поверхности пятна. При этом были опробованы две различные модели испарения, одна из которых предполагает, что испарение нефти регулируется в основном приграничным воздушным слоем, а вторая исходит из предположения, что испарение нефти носит в большей степени диффузионный характер. Расчёты показали, что выбор модели испарения существенно влияет на изменения толщины пятна и суммарного объёма нефти с течением времени.
Ключевые слова: квазилинейное уравнение конвекции-диффузии, метод конечных разностей, генерация разностных схем, базис Грёбнера, испарение нефти.
1. Введение
Для моделирования эволюции нефтяного пятна на поверхности водоёма обычно используется квазилинейное уравнение конвекции-диффузии
дН- -V—- и—+ (—(Б К2—^ + —(Б Н2—м (1)
дЬ дх ду \дх \ дх / ду \ ду)) ' ' ' '
где Н = Н(х,у, — толщина пятна в момент £ в точке (х, у), V и у — скорости дрейфа пятна, Б = д(рш - рац)/(раи • с/), д, Cf — константы, рш и р0ц — плотности воды и нефтепродукта, — функция, позволяющая учитывать наличие процессов деструкции нефти (испарение, эмульгирование и пр.) [1]. Здесь будет учитываться только испарение, поэтому полагаем / = -кеу(¿)Н, где кеу(¿) — скорость испарения [2]. Используя метод расщепления по пространственным переменным, приходим к серии одномерных задач вида
дН дН д / 12дН.
дн = -дхх + дГх(он2ш) - (2)
Статья поступила в редакцию 30 декабря 2013 г.
Авторы выражают благодарность профессору В. П. Гердту за ценные советы и поддержку. Вклад одного из авторов (Ф.Б.Б.) был частично поддержан грантом 13-01-00668 Российского фонда фундаментальных исследований.
Дополнив уравнение (2) граничными и начальными условиями, решение будем искать методом конечных разностей. Введём в [0,1] х [0,Т] равномерную сетку. Зададим шаг по времени Дt, по пространству Дх, количество узлов по х обозначим N, по t — Т. Обозначим К™ = К(гДх,пД£), г = 0, N, п = 0, Т, ^Дж = I, ТДt = Т'. Отметим, что обычно для аппроксимации конвективного члена уравнения (2) используются центральные разности или разности против потока [3]. В случае наличия сильного конвективного течения используются аппроксимации 3-го, 4-го порядка [2,4]. Диффузионный член обычно аппроксимируют со вторым порядком по формуле
(^+1/2^+1 - (а^1/2 + а^1/2) К + аК) (Дх)2
D ^-'--~-"> (3)
где а = К2, а'у+1/2 = (а-+1 + а% )/2, а1?_1/2 = (а% + а-г_1)/2, а = К2 [2,5]. В случае,
когда к = п + 1, к1 = п +1 (п +1 — новый временной слой), для осуществления расчётов необходимо пользоваться итерационными методами.
2. Генерация разностных схем методами компьютерной
алгебры
Применим к (2) метод автоматической генерации разностных схем, предложенный в [6] и реализованный в виде пакета LDA для Maple [7]. Для этого преобразуем (2) к виду
dh_ dh Did2^
~dt = — dX + У ' ( )
где ш = h3, f = 0 (испарения нет). Согласно [6] преобразуем (4) к интегральному закону сохранения. В зависимости от применяемых формул численного интегрирования меняется конфигурация узлов по х, получаются явные схемы (при нижних прямоугольниках), неявные (при верхних прямоугольниках), схемы типа Кранка-Николсона (при формуле трапеций). Например, можно получить неявную схему вида:
h«+i - ц h^1 - h-j1 + D (h^1)3 - 2(h?+1)3 + (h?_+1)3 (5)
At Ах 3 (Дх)
при другом же способе интегрирования
h+ + h-+1 - (h-+1 + h?) _ -v h+ - h?-1 + D (h-+1)3 - 2(h?)3 + (h--1)3
Дt 2Дх 3 (Дх)2
Замечание. В схеме (5) имеется нелинейность на неизвестном слое, поэтому необходимо пользоваться итерационным методом. Для расчётов в данной работе использовался метод Ньютона. В качестве начального приближения выбиралось решение на предыдущем слое.
При увеличении количества точек, задействованных в формулах численного интегрирования получаются более сложные, а также многослойные схемы, например схема Ричардсона [3] для квазилинейного случая.
Также согласно [6] для генерации разностных схем мы можем составить входную для ЬБЛ систему линейных разностных полиномов так: первый полином получаем, дискретизуя (4) на сетке; остальные полиномы получаем, численно интегрируя соотношения / Кг&Ъ = К(х, ¿2) — К(ж, t1), f шх&х = ш(х2, ¿) — ш(х1, ¿),
/ шххёх = шх(х2, ^ — шх(х1, ■£). Положим х2 = (г+ 2) • А, ¿2 = (п + 1) • А1, для интегрирования / Нг&Ъ используем правило Лакса (как в [6]). Комбинируя для / шх ёх, / шххёх формулу средних и трапеций, можно получить 8 различных схем (3 из них одинаковы). Например,
2НГ-+21 - Н™+1 - Н"+3 _ _^ _ Н"+3 - Н"+1 + Б'п ШГ+4 - 2Ч™+2 +
2Аг 2Ах (2Ах)2
2(Н+ + Н++21) - (Н+1
+ К+3 + Н"+5)
4Аг
=
К+4 — К+2 + Б Ч"+6 + — (Ч™+4 + Ш1+2)
2А х
(2Ах)3
Анализ полученных данным способом шести схем с помощью пакета РБЛ [8] показал, что все они имеют порядок аппроксимации 0(АЬ, (Ах)2), а также обладают одинаковой схемной вязкостью, определяемой членом ((Ах)2 • НХХ)/2АЬ. С помощью пакета ЬБЛ можно получать различные разностные схемы и для уравнения (1) непосредственно. Введя промежуточный временной слой и составив систему линейных разностных полиномов по аналогии с тем, как это делалось в [6] для получения двухшаговой схемы Лакса-Вендроффа, можно получить схему расщепления.
3. Численное сравнение конечно-разностных методов
Согласно [5] уравнение (2) при V = 0 имеет точное решение вида
Н аиЬо^-!
(2 • с\/с • £ — х)/Б, если х < с • £
0,
если х > с - Ь
(6)
где Н(х, 0) = 0, Н(0, г) = Н0 • уД, Н0 = ^/2 • с2/Б, с = со^ > 0.
Ь(м)
15.215.11514.914.814.714.614.514.414.3-
О 20 40 60 80 100 120 140 X М
_ Численное решение, полученное с помощью кубической аппроксимации диффузионного члена Автомодельное решение_
100 150
I (мин)
Рис. 1. (а) Численное решение по схеме (5) и точное решение (6) при £ = 2Дt
(б) Относительная погрешность
Исключим из расчётной области точки разрыва —, выбрав с так что с-Т ■
ох
Дг > N ■ Дх, а именно положим с = 600, Т = 300, Д£ = 1 (мин), N = 300, Дх = 1 (м), р0ц = 827 (кг/м3). Приведём расчёты по неявной схеме (5) при V = 0 (рис. 1). А также, заменив в (5) кубическую аппроксимацию диффузионного члена на стандартную (3) при к' = п, проведём расчёты при тех же параметрах (рис. 2). Полученные графики показывают, что обе схемы хорошо аппроксимируют точное решение на том участке, где оно непрерывно. Относительная погрешность Е =|| К — Иаию || / || Каию || быстро убывает с увеличением ¿, причём при использовании (5) она убывает быстрее.
Ь(м)
_ Численное решение, полученное с помощью стандартной аппроксимации диффузионного члена Автомодельное решение_
100 150 200 t (мин)
Рис. 2. (а) Численное решение, полученное при использовании (3), и точное решение (6) при £ = 2Д£. (б) Относительная погрешность
Замечание 1 (об устойчивости). На основе численного эксперимента можно сделать вывод, что (5) устойчива. К такому же выводу можно прийти, применив к (5) метод замороженных коэффициентов. Численные эксперименты с явной схемой, аналогичной (5), показали, что при условиях, указанных в [5] для случая стандартной аппроксимации диффузионного члена, она устойчива.
4. Сравнение различных моделей испарения нефтяного
пятна
Осуществим расчёт по (5) изменений толщины нефтяного пятна в одномерном случае в предположении, что присутствует испарение. Для этого добавим к (5) член к"г+1 ■ К7+1 . Для пересчёта р0ц^) будем пользоваться формулами из [4]. Для расчёта кеу применим 2 методики: 1) из [9], основанную на предположении что ветер существенно влияет на испарение нефти; 2) из [1], отрицающую это предположение на основе проведённой автором [1] серии экспериментов. Результаты расчётов представлены на рис. 3.
5. Заключение
Для численного решения уравнений (1), (2) было получено несколько разностных схем с помощью пакета ЬБЛ. Сравнение с автомодельным решением показало лучшую точность схемы (5) по сравнению с стандартной схемой без итераций.
0.9
V (тл3)
0.6
-Испарения нет 4 Формулы 2й гуппы
° Формулы 1й группы_
h(m)
0.0005-
■ 4|4itifitit|titif И" I I | I I I I | I I I I |
16500 17000 17500 18000 18500 19000 19500
-Испарения нет ■ Формулы 1й гуппы
♦ Формулы 2й группы_
Рис. 3. а) Изменения суммарного объёма б) Изменения толщины
Расчёты толщины нефти с учётом испарения показали, что выбор модели испарения 1) или 2) существенно влияет на изменения толщины пятна и суммарного объёма нефти с течением времени.
Литература
1. Fingas M. Oil Spills Science and Technology: Prevention, Response and Cleanup. — Burlington: Elseivier, 2011.
2. Tkalich P. A CFD Solution of Oil Spill Problems // Environ. Modeling and Software. — 2006. — Vol. 21. — Pp. 271-282.
3. Роуч П. Вычислительная гидродинамика. — М.: Мир, 1980. [Roach P. Computational Hydrodynamics. — Moscow: Mir, 1980. — (in russian). ]
4. Zadeh E., Hejazi K. Eulerian Oil Spills Model using Finite-Volume Method with Moving Boundary and Wet-Dry Fronts // Model. Simul. Eng. — 2012. — Vol. 2012. — P. 33.
5. Тихонов А. Н., Самарский А. А. Уравнения математической физики. — М.: Наука, 1977. [Tikhonov A. N., Samarskii A. A. Equations of Mathematical Physics. — Moscow: Nauka, 1977. — (in russian). ]
6. Gerdt V. P., Blinkov Y. A., Mozzhilkin V. V. Grobner Bases and Generation of Difference Schemes for Partial Differential Equations // Symmetry, Integrability and Geometry: Methods and Applications. — 2006. — Vol. 2. — P. 26.
7. Gerdt V. P., Robertz D. Maple Package for Computing Grobner Bases for Linear Recurrence Relations // Nuclear Instruments and Methods in Physics Research. — 2006. — Vol. 2. — P. 26.
8. Gerdt V. P., Blinkov Y. A. On Computer Algebra-Aided Stability Analysis of Difference Schemes Generated by Means of Groobner Bases // Computer Algebra and Differential Equations, Acta Academiae Aboensis, Ser. B. — 2007. — Vol. 67, No 2. — Pp. 168-177.
9. Stiver W., Mackay D. Evaporation Rate of Spills of Hydrocarbons and Petroleum Mixtures // Computer Algebra and Differential Equations, Acta Academiae Aboen-sis, Ser. B. — 1984. — Vol. 18. — Pp. 834-840.
UDC 531.383:532.516 MSC 65M06, 76D05
Mathematical and Computational Oil Spills Modelling
A. V. Trepacheva, Ph. B. Burtyka
Computer Center of Southern Federal University build. 2, 200/1, Stachki av., Rostov-on-Don, Russia, 344090
The model of oil spill evolution on the water surface based on quasi-linear convection-diffusion equation is considered. To solve the latter the method of finite differences is used. Two-dimensional problem is reduced to one-dimensional by alternating direction method. For one-dimensional case traditional difference approximations are briefly discussed. And also difference schemes obtained by computer algebra based method of automatic generation of difference schemes are considered. For generated schemes the order of approximation and linear numerical dissipation are estimated. For explicit schemes stability conditions are briefly discussed.
In the absence of convection a numerical comparison of traditional implicit difference scheme and one implicit scheme generated automatically with similarity solution of quasilinear diffusion equation is carried out. This comparison shows that generated implicit scheme allows to obtain the relative error less than for traditional scheme.
Using obtained implicit scheme evaluation of oil slick thickness changes in the presence of oil evaporation is done. Two different models were used to estimate oil evaporation rate. The first one is based on hypothesis that oil evaporation is regulated by air boundary layer. The second assumes that oil evaporation is regulated by oil diffusion. Calculations show that the choice of model essentially influences on oil slick thickness and oil total volume changes in time.
Key words and phrases: quasilinear advection-diffusion equation, finite difference method, Grorbner basis, difference schemes generation, oil spills, oil evaporation.
