CC BY
13УДК 511
МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНОЙ ФИЛЬТРАЦИИ, ВЫЗВАННОЙ
РАЗРАБОТКОЙ МЕСТОРОЖДЕНИЙ
С. В. Шешенин1, Э.Р. Какушев2, Н.Б. Артамонова3
В работе рассматривается моделирование осадки поверхности земли в районе нефтедобычи. Используется упругий режим фильтрации. Численное моделирование данного явления осуществляется путем дискретизации краевой задачи по пространственным переменным с помощью метода конечных элементов и разностной дискретизации по времени. Численный алгоритм реализован в виде пакета программ. В демонстрационных целях получены решения задач об откачке жидкости из пятислойного грунта.
Ключевые слова: фильтрация, упругий режим фильтрации, скважина.
In this paper, the subsidence of a ground surface near an oil-producing region is simulated. The elastic mode of filtration is used. A numerical simulation is performed by the discretization of the boundary value problem in spatial variables with the aid of the finite element method and by the finite difference discretization in time. The numerical algorithm in use is implemented as a software package. The problem of pumping a fluid from a five-layer soil is solved as an illustrative example.
Key words: filtration, elastic mode of filtration, well.
Достаточно общей постановкой линейной задачи нестационарной фильтрации является модель Био [1], приводимая к системе дифференциальных уравнений [2]
(CijkiUk,i)j = p,i,
(Kij \ дик,k , д i,j,k,l = 1,2,3. (1)
Уравнения выводятся из закона Дарси, уравнений неразрывности и уравнений равновесия при допущении, что деформации грунта являются малыми [3]. Система дифференциальных уравнений (1) сформулирована относительно вектора перемещения упругого грунта u и изменения давления жидкости p, которые возникают в процессе фильтрации (под изменением давления подразумевается разность между текущим давлением и давлением, которое было в жидкости до начала откачки). В приведенных уравнениях C представляет собой тензор эффективных модулей насыщенного грунта при нулевом давлении жидкости в порах, K — тензор коэффициентов фильтрации, а q — интенсивность источников и стоков. В уравнениях n есть объемная пористость грунта, в — коэффициент сжимаемости жидкости, j — удельный вес жидкости.
При использовании в уравнениях метода конечных элементов (МКЭ) возникает проблема, связанная с необходимостью удовлетворения условия LBB (Ladyzhenskaya-Babuska-Brezzi) [4]. Преодолеть эту проблему можно разными способами, например с помощью более сложных конечных элементов. Второй способ заключается в развязывании уравнений в системе (1). Развязывание можно осуществить на физическом уровне. Этот метод используется в рамках теории упругого режима фильтрации [5-7].
Модель упругого режима выводится с помощью упрощающих допущений в исходной задаче. В случае бесконечной однородной изотропной среды применим к первому уравнению системы (1) оператор дивергенции div [(Л + /) grad div u + /Дu — grad p] = 0. Очевидное преобразование приводит к уравнению Др = 0, где р = div u — ap, а =
(Л + 2/)-1. Если на бесконечности (на практике при достаточном удалении от месторождения) возмущения равны нулю (на практике малы), то для полученного уравнения имеем граничное условие р = 0. Решение уравнения с этим граничным условием единственно и равно
1 Шешенин Сергей Владимирович — доктор физ.-мат. наук, проф. каф. механики композитов мех.-мат. ф-та МГУ, e-mail: [email protected].
2 Какушев Эльдар Рамазанович — асп. каф. механики композитов мех.-мат. ф-та МГУ, e-mail: [email protected].
3 Артамонова Нина Брониславовна — канд. геол.-минерал. наук, ст. науч. сотр. каф. инженерной и экологической геологии геолог. ф-та МГУ, e-mail: [email protected].
нулю. Следовательно, во всей области ё1у и = ар. Подстановка последнего равенства во второе уравнение системы (1) приводит к развязыванию системы. Второе уравнение системы становится уравнением относительно возмущения давления:
(К \ др
сПу ( — • егас! р + а = (а + пв) —. \7 / д
При этом предполагается, что функции д, р и и достаточно гладкие. Это не противоречит тому, что в численном решении участвуют сосредоточенные стоки, поскольку последние используются уже в конечномерной задаче. Эти упрощения приводят к несвязанным уравнениям
(Cijkiuk,i),j = p,',
Kij „ ^ ^ ( 1 o\9p i,j,k,l = 1,2,3. (2)
Ut4,+9 = UT
Теперь рассмотрим случай, когда жидкость течет в одном из пластов, ограниченных сверху и снизу непроницаемыми пластами. Отметим, что такой случай типичен для практики. Поскольку толщина пласта мала по сравнению с его размерами в плане, то можно считать фильтрацию плоской и пренебречь вертикальной компонентой перемещения U3. Это допущение делается только при отыскании поля давления. На том же основании ui, U2 и p зависят только от xi, Х2. Следовательно, первое уравнение системы (2) становится двумерным и к нему можно применить то же обоснование, что и для бесконечной изотропной трехмерной среды. Таким образом, мы приходим к заключению, что для проводящего пласта можно решить отдельно второе уравнение системы (2), сведенное к двумерному. Затем поле перемещений можно определить из решения трехмерной краевой задачи для первого уравнения системы (2) уже при известном давлении жидкости p. Конечно, приведенное обоснование и схема решения сохраняют силу в случае фильтрации в нескольких изолированных проводящих пластах.
Упрощения, связанные с упругим режимом, позволяют решать уравнения системы независимо. Теперь сложность решения линейной алгебраической системы на каждом временном шаге определяется для первого уравнения только разбросом компонент тензора модулей упругости, а для второго — анизотропией компонент тензорного коэффициента проницаемости.
Дискретизация по времени осуществляется заменой производной по времени во втором уравнении
д (•)
системы (2) на разностную производную назад на каждом шаге tm = тт: = ((•)т — (•)m_1)/r-
Для численного решения задачи (2) был разработан пакет программ на языке FORTRAN. С его помощью вычисляются перемещения грунта и изменения порового давления в районе месторождения, которое представляет собой прямоугольный параллелепипед. Оно моделируется слоистой средой с горизонтальным залеганием слоев. В программе количество слоев и скважин может быть произвольным. Скважины могут находиться в любом из слоев. В одних скважинах можно задать откачку жидкости, а в других — закачку. Также можно задать время начала работы каждой скважины. Таким образом, программа может применяться для оценки экологических последствий эксплуатации месторождения.
Приведем в качестве примера численный расчет перемещений грунта и изменения порового давления нефти в районе месторождения. Месторождение моделируется 5-слойной средой. Нумерация слоев ведется снизу вверх. Свойства слоев представлены в таблице. Нефтеносным является третий слой. Протяженность моделируемого месторождения следующая: длина 2 км, ширина 1,5 км, глубина 1 км; коэффициент фильтрации нефтеносного слоя 0,1 м/сут., пористость 0,1, сжимаемость нефти 3 • 10"10 1/Па, удельный вес жидкости 9800 Па/м. Откачка осуществляется из четырех скважин, находящихся в 3-м слое. Скважины полностью пронизывают нефтеносный слой. Интенсивность откачки и координаты скважин следующие: xi = 500 м, yi = 500 м, Qi = 100 м3/сут.; Х2 = 1500 м, y2 = 500 м, Q2 = -110 м3/сут.; Х3 = 500 м, Уз = 1000 м, Q3 = -115 м3/сут.; Х4 = 1500 м, y4 = 1000 м, Q4 = 150 м3/сут. Положительная интенсивность откачки соответствует откачке нефти, отрицательная — закачке воды. Скважина 1 включается в начальный момент времени, скважина 2 — через полгода, скважина 3 — через год и скважина 4 — через полтора года. В данной демонстрационной задаче считается, что свойства закачиваемой воды и откачиваемой нефти приблизительно одинаковы.
Граничные условия на давление задаются только на границе нефтеносного слоя. Скважины находятся на достаточном расстоянии от боковой границы рассматриваемой области, поэтому на всей боковой границе нефтеносного слоя давление не меняется: p = 0. Нефтеносный слой перекрыт сверху и снизу
др
запирающими слоями, ввиду чего на этой части границы отсутствует переток жидкости: —— = 0.
дп
Граничные условия на перемещения задаются во всех слоях. Поскольку боковая граница рассматриваемой области находится на достаточном расстоянии от скважин, то здесь перемещения равны нулю: и = 0. На нижней границе области перемещения также отсутствуют: и = 0. Верхняя граница области свободна от нагрузок: CijklUk,lnj = 0, где п — нормаль к границе. Предполагается идеальный контакт между слоями. На рисунке представлено распределение вертикальных перемещений на поверхности земли через 5 лет после начала откачки. В зонах откачки нефти вертикальные перемещения отрицательны: происходит проседание поверхности земли. В радиусе 100 м от скважины 1 поверхность проседает на 4 см, около скважины 3 максимальное проседание составляет 6 см. Рядом со скважинами, через которые закачивается вода (скважины 2 и 4), вертикальные перемещения положительны и равны 5 см. Происходит выпучивание поверхности.
Номер Толщина Модуль Коэффи-
слоя слоя, м Юнга, циент
МП а Пуассона
5 150 500 0,25
4 250 14 000 0,35
3 50 16 000 0,28
2 300 14 000 0,35
1 250 17 000 0,28
Распределение вертикальных перемещений на поверхности земли через 5 лет
после начала откачки
Программа расчета задачи фильтрации создавалась в сотрудничестве с кафедрой инженерной и экологической геологии геологического факультета МГУ. По ней были проведены расчеты, описывающие проседание земной поверхности и изменение давления для одного из нефтяных месторождений Западной Сибири.
Работа выполнена при поддержке гранта РФФИ № 08-05-00578-a.
СПИСОК ЛИТЕРАТУРЫ
1. Biot M.A. Theory of elasticity and consolidation for a porous anisotropic solid //J. Appl. Phys. 1955. 26, N 2. 182-185.
2. Киселев Ф.Б., Шешенин С.В. Разностная схема для задачи нестационарной фильтрации в слоистых грунтах // Изв. РАН. Механ. твердого тела. 1996. № 4. 129-135.
3. Победря Б.Е. Механика композиционных материалов. М.: Изд-во МГУ, 1984.
4. Brenner S.C., Scott L.R. The mathematical theory of finite element methods. Baton Rouge (USA): Springer, 2002.
5. Терцаги К. Теория механики грунтов. М.: Госстройиздат, 1961.
6. Jacob C.E. Flow in groundwater // Engineering hydraulics / Ed. by H. Rouse. N.Y.: Wiley, 1950. 321-386.
7. Gambolati G. Numerical models in land subsidence control // Comput. Methods Appl. Mech. and Eng. 1975. 5, N 2. 227-237.
Поступила в редакцию 07.06.2010
