Управление динамикой растекания нефтяных пятен Текст научной статьи по специальности «Математика»

УДК 551.46+532.596

УПРАВЛЕНИЕ ДИНАМИКОИ РАСТЕКАНИЯ НЕФТЯНЫХ ПЯТЕН

© 2009 г. Н.Н. Загриценко1, Э.Н. Потетюнко'

2

1Южно-Российский региональный центр информатизации Южного федерального университета, пр. Стачки, 200/1, корп. 2, г. Ростов н/Д, 344090, nnz@sfedu.ru

2Южный федеральный университет, ул. Мильчакова, 8а, г. Ростов н/Д, 344090, dnjme@math. sfedu.ru

1Southern Russian Regional Computer Center of Southern Federal University, Stachki Ave, 200/1, build. 2, Rostov-on-Don, 344090, nnz@sfedu.ru

2Southern Federal University, Milchakov St., 8a, Rostov-on-Don, 344090, dnjme@math. sfedu.ru

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

Ключевые слова: нефтяное пятно, осцилляция, идеальная жидкость, растекание, поверхностное натяжение, математическая модель.

In this paper we study a shape monitoring of oil slicks caused by operational and accidental discharges of ship tanks. We assume that spreading of oil slick is axially symmetric and describe its motion by system of the equations for the three-dimential, axisymmetric, non-viscous fluid with the horizontal planar upper and low boundaries. The differential equations for the oil slick radius and thickness functions of time are found and their solutions are obtained by the numerical methods.

Keywords: oil slick, oscillation, non-viscous fluid, spreading, mathematical models.

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

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

Постановка задачи

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

Решение проведем в цилиндрических координатах г, в, г, а растекание нефтяной пленки будем считать осесимметричным У0 - 0. Г У /с О = 0 где V = . V-

и потенциальным:

Qr =

1

dvz д<уд

дв

д<г(

в.

дг

dVr дв

= 0, О, =-

дК ciV,

dz дг

= О,

= 0.

Уравнение неразрывности в цилиндрической системе координат при осевой симметрии имеет вид

дГуЛ diVz

дг дг

Для нашей задачи существует функция тока удовлетворяющая уравнению

д2ц/ д2ц/ 1 ду/ _ Q

dz2 дг2 г дг Граничные условия для (1) зададим следующим образом:

„2

(1)

Vz=-^ = 0,z = 0

г дг

vz = -l-^=h,z=m

г or

Vr=-^- = R(t,z), r=R(t,z) Г dz

Vr=-l-^ = 0r=0

(2)

(3)

(4)

(5)

Здесь R(t), И(() - радиус и толщина нефтяного пятна.

Условие (2) означает, что нефтяное пятно растекается по недеформированной поверхности нижней жидкости, (3) является кинематическим и показывает,

г

1

1

г

г oz

что вертикальная компонента скорости частиц нефти на верхней границе пятна совпадает со скоростью опускания верхней границы пятна, (4) - кинематическое и показывает, что радиальная компонента скорости частиц нефти на боковой поверхности пятна совпадает с радиальной скоростью растекания пятна. Условие (5) отражает условие осевой симметрии.

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

Решая краевую задачу (1) - (5), находим функцию г2г •

тока у/ =--¡1(1), а из условия постоянства объема

2 к

нефтяного пятна в любой момент времени выводим зависимость между R(t) и Иф:

h(t) =

hpRp

A

R2(t) R2(t)?i R

2

(6)

9 4/г(0 4 ' 2/г(0

Потенциал скоростей здесь определяется с точностью до произвольной функции времени / (/).

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

Воспользуемся интегралом Коши [4]

2

'2

+ gz

+ P=F(t) .

Здесь р - плотность нефти; g - ускорение свободного падения; P - давление в жидкости (нефти); V -скорость частиц в жидкости.

Запишем интеграл Коши для точек, лежащих на верхней границе пятна z = h ( . на которой задано внешнее атмосферное давление Р = Р0 = const, и для точек, лежащих на боковой поверхности пятна, на которой, в силу наличия поверхностного натяжения,

задан перепад давлений: Р~-Р+= a(l/Rl+yR2). Здесь I' - давление внутри пятна; I' - внешнее атмосферное давление; Р+ = Р0 = const; а - коэффициент поверхностного натяжения [5]; Rj, -главные радиусы кривизны боковой поверхности. Согласно [4], полагаем FO= 0 .

Уравнения на боковой поверхности и верхней границе выполним в интегральном виде, т.е. имеем 2 жкС

1 Пр

о о

2я-Д<

- Pjdzdq) = О

J \ | Р

О о

+ gz

+ Р Vdrdcp = 0 , z = И ( .

Приравняем эти выражения:

2

'Ж

I

о

i f /dt о V

hi

dr

--K-

-R<

n V

д<р/ v2

P

(7)

=RC

Здесь /I = - объем нефтяного пятна, кото-

рый считаем постоянным (пренебрегая испарением и растворением); А/я = а(]: /?„ - начальная толщина

пленки; Ro - начальный радиус нефтяного пятна; Щ) -скорость опускания верхней границы нефтяного пятна; К (/) - скорость растекания (изменение радиуса) нефтяного пятна.

Зная функцию тока ц/ в изменяющейся по времени области и пользуясь дифференциальной связью между функцией тока и потенциалом скоростей, определяем потенциал скоростей <р в этой области

-А(0 + т77тА(0 + /о(0.

(далее, чтобы не загромождать выкладки, аргумент t будем опускать).

Подставив в (7) выражения для д(р/Б1, V2/2 и

приняв приближенно |__() = Л , 1//^2|г_0 = 0 (в силу

малой толщины пленки), проинтегрируем выражение (7). Получим

R 2 (R

2 { 4

h2h

-h ■

(h}

KhJ

2

-- + Rh ■

(h

■ + g =

(8)

В (8) И и R связаны между собой через объем (6). Подставим выражение для И через Я в уравнение (8). Для этого нам понадобятся производные от И по вре-

2$ • R

мени: h =--R , h = -2а0 —— + 6а0

Я3 'Я3 Я"

Получим дифференциальное уравнение, определяющее закон изменения радиуса нефтяного пятна с течением времени:

RRL

,4 •

r\ rL+^^L-£т_flL R 2

4 3R7 R R4

2aQ 3a0

R

, aQg aQ g , v 2aQ J J aQ r,V

R

R3

R

R

P j

(9)

где у = ajp.

Исследование полученного дифференциального уравнения

Как отмечают в [3], при достижении некоторого радиуса пятно начинает осциллировать. Проверим, будет ли появляться осцилляционный режим в нашей модели, и если да, то каковы условия его появления.

Уравнение (9) представляет собой дифференциальное уравнение 2-го порядка. Для того чтобы исследовать его решения, сведем его к системе дифференциальных уравнений первого порядка заменой

И(1) = л'(7) (далее будем считать: /0 = 0, /'0 = О [4]).

Получим нелинейную автономную систему:

2

z

Г

3

2

2

2

r

R О st

s О

12

8ao-12a0R6-12ao R3 J

R 3R

3j??9 +8<3q -\2a0R6 -\2а^Ъ ^

(10)

J

ку <j n л и особые точки, являющиеся решением

Найдем положения равновесия системы (10), приравнивая правые части уравнений системы к нулю и учитывая, что s C 0. Получим нулевую особую точ-

{,<> 0

\ROo

\s=0

системы < 3 2 •

Чтобы найти эти ненулевые особые точки, необходимо решить кубическое уравнение

з ^

аоgR - 2aoVR - a0g = 0 •

Введем безразмерный радиус пятна соотношением

R = RqR и преобразуем получившееся кубическое

Р У V 2 р

Так как радиус не может быть отрицательным, нам подходит только 2-й корень, т.е. искомой особой точкой системы (10) будет

5 = = 0

R = ro= Ro

2 а

pgR0

pgRoh0

4 а

\

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

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

дующим образом:

. ¡i = a$+brj

i = R->o

i . где 1

\J] + dt] [t] = s-s0

Вычисляя коэффициенты a, b, c, d системы, полу-

чим 1 , . > . где

n

N^yilR3- (0g (0 -R3 }la0]R; AlfiyiR9 +8ag -12a0R6 -12a$R3 .

Очевидно, осцилляционный режим имеет место, если М У 0 и знаки ^'^иМ^ различны.

Сначала определим, при каких начальных данных будет выполняться условие М^^О, затем проанализируем знаки М^ и Ж1' ^ Дальнейший анализ проведем численно, используя пакет Мар1е.

Зададим конкретные начальные значения Я() =10.

Ло = 1, фиксируя таким образом значение £ = ОД. Найдем те значения р, при которых знаменатель обращается в 0. Корни сначала локализуем графически (рис. 1), а затем в окрестности точек пересечения с осью координат используем численные методы.

М<о1 0.4-

0,3 ß

Рис. 1

уравнение к виду Я 3 - {а~/ pgRо2 (о/я0 У 0 .

Обозначим коэффициент при R буквой Д свободный член - е (малый параметр). При фиксированных Р и е «0 решение строим в виде рядов по степеням е. Я = + + + +■■■■

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

Решение будем искать на интервале 0 < р < 1, так как из физических соображений р> 0, а при Р> 1 корней не будет (как видно из рис. 1). Получим, что при £ = 0,1 знаменатель имеет 2 корня: Д =0,11285837)2 и Р2 = 0,4440629944. При 0< Р <0,11285837)2 и Д > 0,4440629944 знаменатель >0. Далее таким же образом определим корни числителя. Числитель будет всюду отрицателен, следовательно, искомые условия осцилляции при е = 0,1: 0< Р <0,11285837)2 и Д > 0,4440629944. Далее построим фазовые траектории для фиксированных значений £ и р (£ = 0,1, Р = 0,05 ) (рис. 2).

Рис. 2

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

ляция, воспользовавшись выражением

2

(плотность нефти, ускорение свободного падения и

начальный радиус нефтяного пятна - величины известные).

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

нефти

а = ßgRо Р 2

,2

можно варьировать за счет обра-

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

Литература

1. Fay J.A. The spread of oil slicks on a calm sea // Oil on the sea. N.Y., 1969. P. 53-63.

2. Fay J.A. Physical processes in the spread of oil on a water surface // Proc. of Joint Conf. on prevention and control of oil spills. Washington, 1971. № 8. P. 1-23.

3. Аникеев В.В., Мишуков В. Ф. Формирование и разрушение пятна нефти на поверхности моря после аварийного разлива // Комплексные исследования проблем антропогенного загрязнения океана. Владивосток, 1981. С. 108-116.

4. Кочин Н.Е., Кибель И.А., Розе Н.В. Теоретическая гидромеханика. Ч. 1. М., 1963. 584 с.

5. Левич В.Г. Физико-химическая гидродинамика. М., 1959. 700 с.

Поступила в редакцию

29 декабря 2008 г.