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

2015 ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА Сер. 1. Том 2(60). Вып. 4

МЕХАНИКА

УДК 532.7

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

А. Н. Бестужева1, А. Л. Смирнов2

1 Петербургский государственный университет путей сообщения, Российская Федерация, 190031, Санкт-Петербург, Московский пр., 9

2 Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7-9

Рассматриваются двумерные и трехмерные задачи о распространении диффундирующего вещества на водной поверхности и в толще воды. В статье предложено аналитическое решение краевых задач для уравнения диффузии в неограниченных областях при начальном условии специального вида. Проанализирована область концентрации диффундирующего вещества выше «порогового». Рассмотрены случаи распространения диффундирующего вещества по свободной поверхности и на дне водоема. Аналитические решения задач получены с помощью метода Фурье с последующим разложением произвольной функции по функциям Бесселя и полиномам Лежандра. Построенные аналитические решения сравниваются с численными решениями краевой задачи, полученными в пакете Ма^ета^са. Исследована зависимость размера «пятна» загрязнения от времени, а также влияние геометрических и физических параметров на величину радиуса «пятна». Рассмотренные математические модели имеют важное прикладное значение в проблеме защиты окружающей среды при возникновении аварийных ситуаций на морских судах. Библиогр. 12 назв. Ил. 5.

Ключевые слова: диффундирующее вещество, уравнение диффузии, пятно загрязнения.

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

Сходные вопросы были рассмотрены ранее в работах, посвященных решению уравнения теплопроводности в цилиндрических и сферических координатах с граничными условиями первого, второго и смешанного рода [4, 5]. В статье [6] решались одномерное и двумерное уравнения диффузии для стержня и круга с различного ви-

* Работа выполнена при финансовой поддержке СПбГУ (НИР 6.15.820.2013) и РФФИ (гранты 13-01-00523 А и 15-01-06311 А). Часть результатов исследования была включена в доклад на конференции Седьмые Поляховские чтения (Санкт-Петербург).

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

В работе [7] получены уравнение диффузии как в общем виде, так и для изотропной среды в цилиндрических телах, а также выражение для коэффициента диффузии, исходя из физических соображений. Методы решения задач для уравнения диффузии в случае трехмерной ограниченной области подробно описаны в [5].

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

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

Описание распространения диффундирующего вещества сводится к решению краевой задачи для уравнения турбулентной диффузии [2, 7]. Это уравнение не учитывает конвективный перенос диффундирующего вещества, но должно учитывать как адвекцию (движение жидкости относительно выбранной системы координат), так и турбулентный тип движения жидкости. Для несжимаемой жидкости в отсутствие источников (стоков) диффундирующей примеси уравнение турбулентной диффузии для функции концентрации диффундирующего вещества с = с(х, у, г,Ь) принимает вид [2]

дс

— + иУс = Б Ас + У{КУс),

где Д = + Щ2 + -§¿2 — оператор Лапласа, В — коэффициент (молекулярной) диффузии, К — тензор коэффициентов турбулентной диффузии, и — вектор скорости суммарных поверхностных и ветровых течений.

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

дс „ (д2с д2с\ д дс

где коэффициент К2 может зависеть от вертикальной координаты г и времени £ [2]. Это уравнение и его модернизированные и упрощенные варианты широко используются для решения диффузионных задач (например, [1, 2, 8]).

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

3. Распространение диффундирующего вещества по поверхности воды. Распространение диффундирующего вещества в двумерном приближении при

условии, что влияние течения (поверхностного и/или ветрового) не учитывается, то есть при и = 0 описывается уравнением диффузии, в котором неизвестной функцией является концентрация диффундирующего вещества с = с(х, у, £):

с начальным условием

и граничными условиями

дс (д2с д2с

дЬ \дх2 ду2

с(х,у, 0) = /(х,у)

с(г, ¿) ^ 0 при г ^ то

с(0,г) < то.

(1) (2) (3)

Здесь г = \/х2 + у2, а К = К— горизонтальный коэффициент турбулентной диффузии, который далее считается постоянным.

Предполагаем, что диффундирующее вещество, плотность которого меньше плотности воды, в начальный момент времени заполняет круг с радиусом г = I. Перейдем к полярной системе координат в плоскости, совмещенной со свободной поверхностью, при этом учтем, что в силу изотропии ^ = 0. Запишем уравнение (1), начальные (2) и граничные (3) условия в безразмерных переменных: с1 = г1 = у, Ь1 = -уг, где I и со —характерные значения длины и концентрации, в качестве которых в дальнейшем принимаются радиус «пятна загрязнения» и концентрация диффунди-руюшего вещества в начальный момент. Индекс 1 далее опускаем.

Краевая задача примет вид уравнения

дс д 2с 1 дс д^ дг2 г дг

(4)

с граничными условиями

с(г, ¿) ^ 0 при г ^ то

с(0,г) < то

(5)

и

и

и начальным условием

с(г, 0) = / (г). (6)

Используя метод Фурье, получим общее решение уравнения (4) в виде

сю

с(г,*) = У (ал7о(Лг)+ ВаУо(Аг)) е-хН3,\, (7)

о

где Jo(Лг), У0(Лг) —функции Бесселя первого и второго рода соответственно [9]. Из второго граничного условия (5) получаем Вл = 0, первое граничное условие (5) выполняется автоматически.

Произвольную константу Ал определим из начального условия (6). Подставим это условие в решение (7) и используем интегральную формулу Ганкеля [9]:

с с

/(г) = | Jo(Лг)ЛdлJ /(£ША£)^£,

0 0

откуда

сМ) = У в-л2 . (Лг)ЛЙ^ /

о о

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

при котором

с(г, 0)

1, 0 < г < 1, 0, г > 1,

СЮ 1

'' ,-л2г ''

с(г,£) = е-л .0(Лг)ЛЙЛ / .0(Л£)£С

Поскольку

1

получим функцию концентрации в виде

с(г,£)^У .0(Лг). (Л)е-л2^Л.

(9)

Согласно [9] решение (9) может быть представлено в виде рядов

1

с(гЛ) = — .

V > > и ^

^ (-1)

т / _2

— ] 2^1 ( -то,-то, 2,-^-

(10)

с(г, £)

1

(-1)

4£ т=о (т + 1)!

2^1 (—т, —1 — т, 1, г2)

(11)

Здесь 2^1 (а, 6, с, г) — гипергеометрическая функция.

Заметим, что в [10] приведена иная форма решения задачи с начальным условием вида (8):

= 1/И*/.

(12)

где /0(ж) — модифицированная функция Бесселя.

Сравнение численного решения исходной задачи (4)-(6), полученного с помощью сеточных методов, с аналитическими решениями, полученными по формулам (9)-(12), выявило их полное совпадение. С вычислительной точки зрения формулы, содержащие ряды, более удобны, однако при малых значениях £ и больших значениях г происходит потеря точности при суммировании знакопеременного ряда больших величин. Формула (9) не содержит особенности при £ = 0, но предполагает интегрирование по бесконечному интервалу.

На рис. 1 построена зависимость концентрации диффундирующего вещества с(г, £) от времени и расстояния.

1

т

Рис. 1. Зависимость концентрации диффундирующего вещества с(г, 4) от времени и расстояния.

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

Построим зависимость радиуса такого «пятна» от времени т(£; ст). Для этого следует найти г из уравнения

с(т,г) = ст,

где с(г,¿) вычисляется по любой из формул (9)-(12). На рис. 2 приведены результаты численного решения этого уравнения в пакете МаШетайса для разных значений ст

(0 < ст < 1).

Рис. 2. Зависимость радиуса «пятна» загрязнения от времени для различных значений ст: 1 — 0.1, 2 — 0.2, 3 — 0.3, 4 — 0.4.

Для ст < 1/2 размер пятна загрязнения сначала растет, достигая максимума, а затем убывает, при ст ^ 1/2 пятно загрязнения монотонно уменьшается. Максимальный радиус пятна загрязнения, определяющийся равенством тт(ст) = г(Ьт; ст), где ¿т —корень уравнения г'(¿; ст) = 0, следует искать численно.

С течением времени концентрация диффундирующего вещества уменьшается. Определим момент времени, когда пятно загрязнения с концентрацией выше, чем ст, исчезает. Для этого следует найти время Т из уравнения с(0, Т) = ст. Используя формулу (12), получим

-41'Г'^

откуда

2

= i _ e-1/(4T)

Т- 41п(1 — ст)' (13)

В реальных приложениях ст очень мало (порядка 10-6 — 10-8). Оставляя в формуле (13) только первый член в разложении по малому параметру ст, получим

Г = г-

Определим время жизни пятна с концентрацией не меньше, чем ст, и радиуса не меньше, чем г*, т.е. найдем Т из уравнения с(г*,Т) = ст или

1 ¡Л £ -аг*)2+£2)Я4ти

2./о Т~ 0

Поскольку нас интересуют малые значения cm и, как следствие, большие значения T, разложим подынтегральную функцию в ряд по малому параметру ^ = 1/T, оставив два первых члена, и проинтегрируем. Получим

2

Из этого уравнения находится величина ^ и, следовательно, время жизни пятна T. Учет последнего члена имеет смысл, если он по порядку близок к первому, то есть для пятен с большим радиусом (г*)2^ = O(1).

Рассмотрим пример, в котором толщина пятна загрязнения является постоянной и равной 10-4 м. Полагаем значения параметров l = 0.1 м, K = 102 м2/с, cm = 10-6 кг/м3, со = 680 кг/м3, что соответствует значениям поверхностных концентраций cJm = 10-10 кг/м2 и с0 = 0,068 кг/м2. При этом время существования пятна загрязнения составляет T « 5 часов.

4. Распространение диффундирующего вещества в толще воды. Осе-симметричный случай. Теперь рассмотрим следующую задачу. Пусть диффундирующее вещество с плотностью выше, чем плотность воды, находится в начальный момент на плоском дне водоема. Введем систему координат XOY на дне водоема, совместив начало координат с центром круга, являющегося проекцией области, которую занимает диффундирующее вещество в начальный момент времени. Ось Z направим вертикально вверх. Процесс диффузии может быть описан уравнением

дс т, д2с д2с д2с

9t=K'W+K>W+K'd*- (14)

с

и

Будем считать, что водная среда трансверсально изотропна, то есть К = Кх = Ку и Кг ^ К. Проведем масштабирование переменных: х = х\, у = ух, г = Кг\. Опуская индекс 1, получим уравнение, которое после перехода к сферическим координатам примет следующий вид:

дс

1 д

dt

r2 дт

дс

дт

1 1 д sin в т2 дв

■ адс 811119дв

+

1 д2с

sin2 в дф

2

Учитывая изотропию среды в горизонтальной плоскости, рассмотрим осесимметрич-ную задачу, полагая -щ = 0. Тогда краевая задача будет состоять из уравнения диффузии

д2с 2 дс 1 д2с 1 со ев дс

дс dt

дт2 т дт т2 дв2 т2 sin в дв'

граничных условии

с(т,в,г) ^ о с(0, в, t) < то, дс

при т ^ то,

при

0 = 1

и начального условия

с(т, в, 0) = f (т, в).

Разделяя переменные в виде c(r, в, t) = T(t)R(r)©(e) и используя в качестве констант разделения Л и ц — волновые числа, приходим к системе обыкновенных дифференциальных уравнений [3]. Следуя [3], положим ц = m(m + 1). Используя [9, 11], получим решения этой системы:

T (t) = C1e-x2t,

©(в) = C2Pm(0+ CsQm(^),

R(r) = CAjm(Xr) + C5 ут(Лг),

Где pm ( x) и Qm(x) —функции Лежандра первого и второго рода, jm(x) и ym(x) — сферические функции Бесселя первого и второго рода, а £ = cos в.

Из условия ограниченности решения при в = 0 следует, что C3 = 0. Из условия дР^в ^ = 0 ПРИ ^ = § следует, что то— целое число, причем четное: то = 2п. Из условия ограниченности решения при r = 0 следует, что C5 = 0. Условие c(r, e,t) ^ 0 при r ^ то выполняется автоматически.

Выпишем общее решение краевой задачи:

с(т, в,г) =

то

Е

n=0 '

с+то

Ах,и32и(Хт)Р2п(cosв)е ld\.

(15)

Константы A\,n находятся из начального условия c(r, в, 0) = f (r, в). Рассмотрим случай, когда f (r, в) = f (r). Полагая t = 0, домножая равенство (15) на P2m (cos в) sin в и интегрируя по в от 0 до получим, что правая часть равенства отлична от 0 только при n = m, а левая — только при m = 0. Следовательно,

A\jo(^)d\

sinAr Ax^—d\ = f(r).

0

0

0

Раскладывая функцию g(r) = г/(г) в интеграл Фурье, получим

сю сю сю

g(r) = J —^ sin XrdX = — J J g(£) sin sin ArdA,

00

откуда

0

Для начального условия

справедлива формула

с(г,0, 0)= /(r)H 1, 0 > 1 > 1,

/ 2 (sin А , . .

Ax = J-l — -coS\), (16)

и выражение для функции концентрации принимает вид

ТО 2

. . 2 Г f sin A \ e-xt . Л 1Л

c(r,t) = — / —---eos А —-—smAraA. (17)

пг J V A JA

о

Отметим, что для исходной задачи (в размерном виде) такое начальное условие соответствует сильно сплющенному в вертикальном направлении эллипсоиду с полуосями в отношении 1: у/Kz/K.

В [10] приведена другая форма решения этой задачи для рассматриваемого начального условия:

1

1 t ( (r-A)2 (r + A)2 ,

c(r,t) = -X [ — ¿t ¿X. (18)

Значения функции концентрации можно вычислять по любой из формул (17) и (18). Вычислительные трудности возникают только при очень малых значениях ti « 10-3. Например, при вычислении по формуле (17) в окрестности точки разрыва разлагаемой функции c(r, 0) наблюдается явление Гиббса (рис. 3).

При вычислении концентрации c(r, t) для t > 0.01 в формуле (17) в качестве верхнего предела можно брать небольшое значение Ai. На рис. 4 проводится сравнение значений функции концентрации, вычисленной по формулам (17) и (18) при ti = 0.01 для 0 < r < 5, причем Ai = 18. При больших значаниях Ai графики практически неразличимы.

Вычислим время существования области загрязнения с концентрацией выше, чем cm. Принимая во внимание, что, как и ранее, малые значения предельно допустимой концентрации приводят к долгому существованию области загрязнения, введем малый параметр е2 = 1/t и разложим подынтегральную функцию в формуле (18) в ряд по малому параметру, сохранив два первых члена:

1 (,+л)2 \ е3А2 е5А2,2 л2.

-—А е — - е — = —--(г2 + Л2).

2 rsJVt V / 2yi 8^ ^ ^

Рис. 3. Явление Гиббса.

Рис.4. Функция с(г, ¿х), вычисленная по формулам (17) (синяя линия) и (18) (красная линия).

После интегрирования получим

£3 £5 /г2 1 Ст = " (уУ + 5

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

1

Т :

62/3с'3П1/3 '

Понятно, что учет размера пятна имеет смысл только, если ге = 0(1).

Отметим, что в отличие от [4], где получены асимптотические оценки, связывающие концентрацию и время жизни пятна, здесь в явном виде приведены первые члены разложения для времени жизни пятна.

В качестве примера используем значения параметров из предыдущего раздела, но возьмем другие коэффициент диффузии К = 10-1 м2/с и начальную концентрацию со = 1300 кг/м3. В этом случае получим, что время существования пятна загрязнения составляет Т « 7 часов.

При начальном условии произвольного вида для получения коэффициента Л\ следует использовать формулы разложения для интеграла Фурье—Бесселя [12].

На рис. 5 значения функции концентрации с(г, £) вычислены по формуле (17) для разных значений времени. С ростом времени вклад больших Л уменьшается и требуемая точность вычислений достигается при небольших значениях верхнего предела интеграла.

c

1.0

0.4

0.8

0.6

0.2

0.5

1.0

1.5

2.0

2.5

3.0

r

Рис. 5. Функция с(г, Ь) в моменты времени Ь = 0.0001 (синяя линия), Ь = 0.01 (зеленая линия), Ь = 0.1 (красная линия).

5. Заключение. Получено аналитическое решение краевых задач для уравнения диффузии в неограниченных областях при начальном условии специального вида. Найдена область диффундирующего вещества с концентрацией выше «пороговой».

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

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

Авторы выражают благодарность заместителю генерального директора по научной работе ЗАО «НПФ "АРГОС"» В. В. Лобынцеву за предложенную для исследования задачу и полезные обсуждения.

Литература

1. Марчук Г. И. Математическое моделирование в проблеме окружающей среды. М.: Наука,

2. Озмидов Р. В. Диффузия примесей в океане. Л.: Гидрометеоиздат, 1986.

3. Кошляков Н. С., Глинер Э. Б., Смирнов М. М. Уравнения в частных производных математической физики. М., 1970.

4. Socolofsky S.A., Jirka G.H. Environmental Fluid Mechanics 1: Mixing and Transport Processes in the Environment. Texas A&M University. 2005.

5. Olver P. J. Introduction to Partial Differential Equations. Springer. 2013. P. 635.

6. Caretto L. S. Solution of the Diffusion Equation. California State University, College of Engineering and Computer Science, Mechanical Engineering Department. January 25, 2015. P. 36.

7. Thambynayagam R. K. M. The Diffusion Handbook: Applied Solutions for Engineers. McGraw-Hill. 2011. P. 2048.

8. Астраханцев Г. П., Меншуткин В. В., Петрова Н. А., Руховец Л. А. Моделирование экосистем больших стратифицированных озер. СПб.: Наука, 2003. С. 363.

9. Бейтмен Г., Эрдейи А. Высшие трансцендентные функции. М.: Наука, 1966. Т. 2. С. 296.

10. Полянин А. Д. Справочник по линейным уравнениям математической физики. М.: Физмат-лит, 2001. С. 576.

11. Бейтмен Г., Эрдейи А. Высшие трансцендентные функции. М.: Наука, 1965. Т. 1. С. 296.

1982. С. 320.

12. Левитан Б. М. Разложение по функциям Бесселя в ряды и интегралы Фурье // УМН, 1951. Т. 6. Вып. 2(42). С. 102-143.

Статья поступила в редакцию 26 марта 2015 г. Сведения об авторах

Бестужева Алла Николаевна — кандидат физико-математических наук, доцент; bes_alla@inbox.ru Смирнов Андрей Леонидович — кандидат физико-математических наук, доцент; a_l_smirnov@mail.ru

PROPAGATION DYNAMICS OF DIFFUSIVE POLLUTANTS ON THE WATER SURFACE AND IN THE WATER

Alia N. Bestuzheva1, Andrei L. Smirnov2

1 St. Petersburg State Transport University, Moskovskiy pr., 9, St.Petersburg, 190031, Russian Federation, bes_alla@inbox.ru

2 St.Petersburg State University, Universitetskaya nab., 7-9, St.Petersburg, 199034, Russian Federation; a_l_smirnov@mail.ru

Firstly, the 2D problem of propagation of diffusing pollutant on the water surface is analyzed. Such model may be used, for example, to study the lifetime of the toxic pollutant spot on the water surface. For isotropic medium the boundary value problem for the diffusion equation is considered, the analytical solution of which may be obtained by means of Fourier method with consequent expansion of the arbitrary function in Bessel functions. The found analytical solution is compared with numerical solutions of the boundary value problem obtained with Mathematica software packages. The dependence in time of the pollution spot size is studied and the effect of geometrical and physical parameters on the pollution spot radius is discussed. Also the 3D problem of toxic pollutant propagation set on the flat bottom is examined. The size of the domain, where the concentration of the toxic pollutant is higher than the maximum permissible concentration, and the dynamics of this domain are studied. Refs 12. Figs 5. Keywords: diffusing pollutant, diffusion equation, toxic pollutant spot.

References

1. Marchuk G.I., Mathematical Models in Environmental Problems (North Holland, Amsterdam, 1986).

2. Ozmidov R. V., Diffusion of contaminants in the ocean (Kluwer Academic Publishers, 1990).

3. Koshlyakov N. S., Smirnov M.M., Gliner E. B., Differential Equations of Mathematical Physics (North-Holland, Amsterdam; Interscience (Wiley), New York, 1964).

4. Socolofsky S. A., Jirka G.H., Environmental Fluid Mechanics 1: Mixing and Transport Processes in the Environment (Texas A&M University, 2005).

5. Olver P. J., Introduction to Partial Differential Equations (Springer, 2013).

6. Caretto L. S., Solution of the Diffusion Equation (California State University, College of Engineering and Computer Science, Mechanical Engineering Department. January 25, 2015).

7. Thambynayagam R. K. M., The Diffusion Handbook: Applied Solutions for Engineers (2011).

8. Astrakhantsev G. P., Menshutkin V. V., Petrova N. A., Rukhovets. L. A., Modelling of ecosystems of large stratified lakes (Nauka, St. Petersburg, 2003) [in Russian].

9. Bateman H., Erdelyi A., Higher transcendental functions 2 (Dover Publications, 2007).

10. Polyanin A. D., Handbook of Linear Partial Differential Equations for Engineers and Scientists (Chapman & Hall/CRC Press, Boca Raton-London, 2002).

11. Bateman H., Erdelyi A., Higher transcendental functions 1, (Dover Publications, 2007).

12. Levitan B.M., "Expansion in Fourier series and integrals with Bessel functions", Uspekhi Matematicheskikh Nauk 62(42), 102-143 (1951) [in Russian].