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

Математическое моделирование загрязнения почвы при разливах нефти

Аварийные и штатные разливы нефти и других жидких углеводородов из различных объектов МТП и ПНГХК наносят значительный ущерб природной среде. В результате постоянных протечек происходит проникновение нефтяных загрязнений в пористую среду почвенного покрова и грунтовые воды, что приводит к деградации почв и растительности и требует значительных финансовых затрат на рекультивацию и другие природоохранные мероприятия.

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

где С - концентрация загрязнения, v - скорость течения по трещинам почвы, D - коэффициент диффузии, X - коэффициент трансформации загрязнений в результате взаимодействия с растительным покровом и в результате химических реакций, z - глубина проникновения жидкого загрязнения в почву.

Упрощенная модель (3.1.1) имеет существенный недостаток, присущий большинству упрощенных моделей: ограниченное число параметров, численные значения которых могут быть получены только в результате многочисленных экспериментов.

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

где v - скорость потока через пористую среду, к - коэффициент фильтрации, Р = p + pgz - суммарный напор, z - высота относительно фиксированного

уровня.

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

где s,s =1-5 - нефтенасыщенность и газонасыщенность пористой среды, р,рх,р,рх - плотности и вязкости соответственно нефти и газа, т - пористость почвы, Р,Р - давления в жидкой и газовой фазе, л - абсолютная проницаемость, k(s),kAs) - относительная проницаемость для жидкой и газовой фазы, g - вектор ускорения свободного падения |g| = 9,8м-с'2.

На границе жидкой и газовой фазы в пористой среде выполняется баланс давлений:

где о - коэффициент межфазного натяжения (между газовой и жидкой фазой), 0 - краевой угол, который жидкая капля составляет со стенкой поры, при этом в условиях полного смачивания 0 = 0, J(s) - безразмерная функция Леверетта.

Для моделирования функции Леверетта используют различные аппроксимации, например формулы (3.1.5) и (3.1.6):

где so=0,l - максимальная нефтенасыщенность пор, Jo =0,034 при (т = 0,03 Пам.

Для моделирования относительной проницаемости также используют многочисленные варианты эмпирических зависимостей, например формулы (3.1.7)-(3.1.9):

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

157

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

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

(3.1.3). При этом в формуле (3.1.4) можно пренебречь давлением в газовой фазе по сравнению с давлением в жидкой фазе. Примем также условие не зависимости от z и 5 численных параметров задачи сг,ц,т,р,п Вместе с начальными и граничными условиями эти упрощения приводят к следующей краевой задаче:

Выполняя дифференцирование, преобразуем (3.1.10) к следующему виду:

В уравнении (3.1.11) будем использовать следующие численные значения параметров: р = 1,210~3 Па-с, <т = 0,03 Па-м, 0 = 60°, я- = 10~'2 м2, /> = 760 кг-м3, g = 9,8 м-с’2, т = 0,2. В результате приведем (3.1.11) к следующему виду:

Для моделирования функции Лаверетта и относительной проницаемости будем использовать упрощенные эмпирические зависимости (3.1.6) и (3.1.9), производные которых имеют следующий вид:

Преобразуем частные производные, входящие в (3.1.12):

В результате краевую задачу (3.1.12) окончательно представим в следующем виде:

Численное решение задачи (3.1.15) будем искать на прямоугольной сетке с размером элементарной ячейки AtxAz. Для аппроксимации производной по времени будем использовать метод Кранка - Николсона:

где Ф^. - правая часть уравнения (3.1.15), к, j - временной и пространственный индексы.

В уравнении (3.1.16) для расчета Ф ,Фы . используется следующая конечно-разностная аппроксимация:

Представим уравнения (3.1.16) и (3.1.17) в форме, необходимой для решения задачи методом прогонки

b (s )Д/

где WJ=sk+lJ, А, = С, = DkJ, BJ=l+2DkJ, D =-L-b—,

2 Az

Tj = + "su)+^r^k>. ~skj)2 + Dkj (j-t +j«-2skJ)

Алгоритм, реализующий метод прогонки, называется алгоритмом Томаса. Рассмотрим основные особенности и сущность алгоритма Томаса, имеющего прямой и обратный ход.

Для решения уравнения (3.1.18) методом прогонки введем два вспомогательных вектора Е и F :

Понижая индекс в (3.1.19) на единицу, получим:

Подставим (3.1.20) в (3.1.18) и разрешим полученное уравнение относительно Wj

Сравнение полученного уравнения с (3.1.19) приводит к следующим рекуррентным формулам для вычисления векторов Е и F :

Расчет векторов Е и F по рекуррентным формулам (3.1.22) называется прямым ходом алгоритма Томаса.

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

Граничное условие первого рода (Дирихле) на левой границе имеет следующий вид:

Перепишем уравнение (3.1.19) для левой границы

Из сопоставления (3.1.23) и (3.1.24) для векторов Е и F получим следующее решение:

Граничное условие второго рода (Неймана) на левой границе имеет следующий вид:

Из сопоставления (3.1.26) и (3.1.24) для векторов Е и F получим следующее решение:

Граничное условие третьего рода (Робина) на левой границе имеет следующий вид:

Из сопоставления (3.1.28) и (3.1.24) для векторов Е и F получим следующее решение:

Из уравнения (3.1.29) следует необходимость выполнения следующего условия: АхфЬг

Рассмотрим теперь граничные условия на правой границе.

Граничное условие первого рода (Дирихле) на правой границе имеет следующий вид:

Граничное условие второго рода (Неймана) на правой границе имеет следующий вид:

Выразим граничное условие для вектора W через вектора Е и F , для чего запишем уравнение (3.1.20) при условии j = N

Приравнивая правые части (3.1.31) и (3.1.32), получаем

Граничное условие третьего рода (Робина) на правой границе имеет следующий вид:

Приравнивая правые части (3.1.34) и (3.1.32), получаем

Как и в предыдущем случае необходимо выполнение условия не равенст-

, Ах

ва нулю знаменателя bN ф-.

Ех, | 1

Приведем программу алгоритма Томаса на m-языке среды MatLab. %TOMAS алгоритм численного решения трехдиагональной % системы линейных алгебраических уравнений методом прогонки

%

function W=Tomas(A,B,C, T,N,L1 ,LN,A1 ,B1 ,AN,BN,dx)

% A,B,C,T - векторы коэффициентов уравнения %o N - количество узлов сетки

%L1, LN - задают вид граничных условий слева и справа % А1, В1 - параметры граничного условия слева % AN, BN - параметры граничного условия справа %dx- пространственный шаг модели

%

E=zeros(l,N); % вспомогательные векторы F=zeros(l,N);

if LI==1 % условия Дирихле слева E(1)=0;

F(1)=A1;

end

if LI==2 % условия Неймана слева E(l)=l;

F(l)=Al*dx;

end

ifLl==3 % условия Робина слева E(l)=-Bl/(dx-Bl);

F(1)=A1 *dx/(dx-Bl); end

% прямой ход i=2:N-l

Den(i)=B(i)-C(i)*E(i-l);

E(i)=A(i)./Den(i);

F(i)=(T(i)+C(i). *F(i-l))./Den(i);

if LN==1 % условия Дирихле справа

W(N)=AN; end

if LN= =2 % условия Неймана справа W(N)=(F (N-l)+AN*dx)/(l -E(N-1)); end

if LN= =3 % условия Робина справа W(N)=(F(N-1) *BN+AN*dx)/(dx+BN-E(N-l) *BN); end

% обратный ход j-N-l:-l:l;

W(j)=E(j). *W(j+l)+F(j); return

Следует обратить внимание на использование оператора "точка", обозначающего в MatLab операцию поэлементного умножения матриц и векторов. Результаты компьютерного моделирования представлены на рис.3.1.1.

Распространение нефтяных загрязнений по глубине почвы

Рис. 3.1.1. Распространение нефтяных загрязнений по глубине почвы: 1-0,1 часа; 2 - 1 час; 3-20 часов; 4-50 часов; 5 - 500 часов

Как следует из графиков (рис.3.1.1) скорость распространения нефти по глубине почвы достаточно мала для выбранных значений параметров почвы. Для моделирования использованы экспериментальные данные по распространению мазута в песках. При абсолютной проницаемости песка я- = 3-10’12 м2 глубина проникновения в один метр достигается мазутом за время около двух часов. Таким образом, результаты компьютерного моделирования свидетельствуют о существенной зависимости скорости распространения нефтепродуктов в почве от состава и влагосодержания грунтов.

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

 
Посмотреть оригинал
< Пред   СОДЕРЖАНИЕ   ОРИГИНАЛ     След >