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

УДК 532.62

Д.А. Масленников, Л.Ю. Катаева, О.Е. Куркина, А.А. Куркин

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

Нижегородский государственный технический университет им. Р.Е. Алексеева

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

Ключевые слова: тонкая плёнка, загрязнения, итерационно-интерполяционный метод, численное моделирование.

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

Задача распространения плёнки может быть сформулирована следующим образом [1]:

—+—(иГ)+ —(уГ) = D

dt dx dy

dx

d2

-(r)V(Г)

+

Го-Г

(1)

В данном уравнении Г - это концентрация ПАВ, Го - начальная концентрация, и и V это проекции внешнего поля скорости, задаваемые в следующей форме: / (х, у, = А соб(^ - £х+ф1 )соб(^ - 1у + ф2 ). Начальные и граничные условия:

dГ

Г(Х,y,t = о)=Го ,Г|,=0 =Го, —

y=У0

= 0, Г(0, y, t ) = rf ^, y, t

(2)

Для учёта периодических граничных условий по оси х решение ищется на области 6л

размером —, то есть периодическое граничное условие заменяется условиями

dT

dx

2л

к ' 5Г

dx

4л

, что хотя и увеличивает время вычислений, но не требует использования

циклическом прогонки.

Реализованный алгоритм состоит из нескольких блоков: Model, Solver2D, SolverlD, Border Condition, TDMA. Model определяет способ дискретизации по времени, вычисляет коэффициенты уравнения и граничных условий для сведения их к обобщённому виду, делегируя решение этого уравнения блоку Solver2D. Solver2D преобразует многомерное уравнение во множество одномерных уравнений гиперболического или параболического вида [2, 3], затем они решаются с помощью SolverlD. SolverlD реализует модифицированную схему итерационно-интерполяционного метода [4]. Он делегирует дискретизацию граничных условий блоку BorderCondition. Итоговая схема (3) - (4), реализуемая данным алгоритмом, имеет вид:

т

x=-

к

к

© Масленников Д.А., Катаева Л.Ю., Куркина О.Е., Куркин А.А., 2011.

At

= D

та2+4—

n + ^ + Г 12

J

n + 1

- i+1, Л/

i/Л Л \ii n "Г n _ j/ n pn )

j-j_(r«i,J + 4rnj + -+,J)J + 6(U+J-+i,J ^-1+^J--1+JJ

(3)

- 2Г"72 +Г"+12 1/ ч

i-1,J ivJ -i+J + f- (Г-1,J + 4-J + Г+1,J - 6Г0)

Ax

2

2т

1 Г(Г" J- + 4-n+1 +

-iv j

Г"J++1 bfa$ + 4Гn+K + -n +■# 11 + 6

-iv j

un гn+У2 -un Гn+У2

ui, J+i- i, J+i ui, J-1+i- i, J-1+i

Ax

(4)

1

Г"+1 -9Г"+1 1 / 1/

_ ^ Ч;-1 2Ч; +Ч;+1 + 1 (г «+/2 + 4Г«+/2 +г «+/2 6Г - V + 2г7г']+1 + г'] + Ч]+1 - 610

где ^ = 0, если п - нечётное, и ^ = 1 в противном случае.

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

(

At < min

1 , 1 5D 5D

Av ^TT^Tj

л

4Ax A '4AyA V Ax2 ' Ay

(5)

где Au, Av -амплитуда скоростей u и v соответственно.

Главное окно программного комплекса содержит в своей верхней части главное меню, в нижней части - строку состояния, а в центральной - MDI зону (Multiple document interface) для работы с окнами построенных графиков. Главное меню служит для управления программным комплексом, а в строке статуса отображается его текущее состояние.

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

Для ввода параметров в программе предусмотрены диалоги. Диалог, представленный на рис. 1, позволяет пользователю задать значения скорости u или v, в зависимости от выбранного пункта подменю «Параметры - Начальные условия». Диалог показывает формулу для скорости с параметрами и поля, содержащие текущие значения этих параметров, которые может изменить пользователь. Параметр к, введённый в этом диалоге для одной из скоростей u или v, автоматически применяется к другой скорости, и определяет размер периода по х.

Рис. 1. Диалог Поля скоростей

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

1

Рис. 2. Диалог Поля скоростей

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

Блок «Параметры отображения» отвечает за визуализацию результатов. Решение показывается в МШ-области. Параметр «Пауза отображения» определяет время, на которое будут остановлены вычисления после каждого отображения. Этот параметр нужно делать ненулевым, если расчёты идут слишком быстро. «Шаг времени для графика» определяет, через какие промежутки времени моделируемого процесса будет производиться отображение. Кроме того, данный параметр влияет на шаг времени в вычислениях, чтобы шаг времени для отображения был кратен шагу время в вычислениях. Параметр «Шагов времени для графика» определяет, сколько раз будет выполнено отображение на графике, прежде чем расчёт будет завершён. После запуска программного комплекса, необходимо задать поля скорости с помощью диалогов, запускающихся из соответствующих пунктов подменю «Параметры -Поля скоростей». Затем с помощью пункта меню «Параметры - Основные» нужно указать размеры расчётной сетки, шаг по времени и параметры вывода результатов на экран. После этого запустить расчёт с помощью пункта меню «Расчёт - Запустить».

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

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

Программный комплекс был протестирован при параметрах, показанных на рис. 1-2, причём скорости и и V равны. На рис. 3 показан начальный профиль концентрации на момент времени 1=\ с.

Рис. 3. Распределение концентрации ПАВ в момент t=1 с

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

Рис. 4. Распределение концентрации ПАВ в момент t=2,5 с

Рис. 5. Распределение концентрации ПАВ в момент t=20 с

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

Библиографический список

1. Пелиновский, E.H. Пленки ПАВ на морской поверхности: препринт Института теплофизики СО АН СССР / E.H. Пелиновский, Т.Н. Талипова. 1990. № 219-90. - 42 с.

2. Яненко, Н.Н. Метод дробных шагов решения многомерных задач математической физики / Н.Н. Яненко. - Новосибирск: Наука, 1967. - 196 с.

3. Катаева, Л.Ю. Особенности дискретизации многомерных нелинейных задач. Наука и техника транспорта. 2008. № 4. С. 13-16.

4. Катаева, Л.Ю. Оценка эффективности итерационно-интерполяционного метода при решении двумерных уравнений гиперболического типа / Л.Ю. Катаева Д.А. Масленников // Актуальные проблемы гуманитарных и естественных наук. 2009. № 4. С. 9-17.

Дата поступления в редакцию 08.02.2011

D.A. Maslennikov, L.Yu. Kataeva, O.E. Kurkina, A.A. Kurkin

MODELING THE MOVEMENT OF A THIN SLICK INDUCED BY PERIODICALLY

EXTERNAL VELOCITY FIELD

The numeric simulation of the movement of the thin slick on the surface of the ocean is considered. Algorithm based on the iteration-interpolation method for modeling dynamics of thin slick induced by external velocity field. . Software system that implements this algorithm was developed.

Key words: thin slick, surface-active substances, ocean, pollution, the iteration-interpolation method, numerical simulation.