CC BY
30
НЕЯВНЫЕ ЧИСЛЕННЫЕ СХЕМЫ В ЗАДАЧАХ ФИЛЬТРАЦИИ В
ПОРИСТЫХ СРЕДАХ
IMPLICIT DIFFERENCE SCHEMES FOR FILTRATION THROUGH POROUS MEDIUM PROBLEMS
2 1 2 C.B. Шешенин , Ф.Б. Киселев , Н.Б. Артамонова
S.V. Sheshenin, F.B. Kiselev, N.B. Artamonova
ФГБОУ ВПО "МГСУ"1, МГУ им. М.В.Ломоносова2
Предложены различные способы дискретизации по времени трехмерной задачи фильтрации. Задача решается в связанной постановке по модели Био. Полученные численные схемы обладают свойством высокой скорости сходимости решения.
Various methods of time discretization for 3d filtration problem are proposed. Coupled equations by Biot's model are solving. Resulting numerical schemes have a good rate of convergence.
Задачи механики грунтов, имеющие приложение к разработке геологических ресурсов, являются востребованными на протяжении не одного десятилетия. Продвижение в этой области видится связанным, в том числе, с решением пространственных задач для областей большого размера и сложной формы. В настоящей работе рассмотрена трехмерная задача фильтрации жидкости через пористый грунт вследствие откачки из скважин. Задача решается в связанной постановке по модели Био [3].
Дискретизация по пространственным переменным осуществляется с помощью МКЭ. Целью статьи является представление различных способов дискретизации по времени. Дело в том, что реализация связанной задачи фильтрации сопряжена с трудностями, вызванными плохой обусловленностью матрицы жесткости. Плохая обусловленность определяется структурой матрицы и большим разбросом модулей слоистого грунта, особенно коэффициента фильтрации жидкости.
Преодоление трудности моделирования фильтрации осуществляется на физическом или математическом уровне. К первому относится теория упругого режима фильтрации. Она основана на «развязывании» уравнений с помощью предположений физического характера. Такой подход использовали К. Терцаги [1], С. Джекоб [5], Д. Гамболати [4]. Этот подход приводит к потере свойства связанности и не позволяет получить, например, эффект Мандела-Крайера. Другой подход заключается в «развязывание» системы уравнений с помощью внутренних итераций на каждом шаге по времени. Такой подход предложен в [2] и применялся в [6] применительно к осесим-метричной задаче откачки жидкости из скважины. Аналогичный метод позднее рассмотрен А. Сеттари и др. [7, 8].
1. Запишем уравнения модели Био, которые связывают перемещения и упругого скелета грунта и давление р в жидкости:
6/2П11 ВЕСТНИК _6/2011_МГСУ
id.lv(С: gradй') — gradp dlv (^-Щ • gradp^ + в — dlvu + пРр, х£К, в — -у
где С - тензор эффективной жесткости насыщенного грунта при нулевом давлении жидкости в порах, к - тензорный коэффициент фильтрации жидкости, р?д - вес единицы объема жидкости, у - плотность источника (стока) жидкости, п - объемная пористость, р - коэффициент сжимаемости жидкости. Точкой в системе уравнений обозначена частная производная по времени. Граничные условия записываются в виде:
о-п = (сик1икЛ -р50)пД х£ I? Й = и0, х е
(2)
, х t Ь" р = Р0, хе I?
• gradp) -n = W0, хЕ S
вместе с нулевыми начальными условиями: и = 0, р = О при t — 0, поскольку задача решается в приращениях относительно состояния статического равновесия. Далее через тензор к будем обозначать отношение к(х~)/р?д.
2. Дискретизация задачи по пространственным переменным основана на вариационных уравнениях, соответствующих системе (1) - (2).
!/т/ gradй ■■ С ■■ gradw dV + /т/ gradp • и/ dV — О
г г А Г (3)
)у gradp ■ к ■ gradq dV + + п$р )qdV — J GqdZG
Первый способ дискретизации по времени основан на замене производных по времени во втором уравнении (3) на разностную производную «назад» на каждом шаге по времени £т = хт:
д(-) (-)т - (•)т~1 ~дГ ~ ~
Применяя к (3) метод Галеркина и используя три-линейные функции формы:
и? = К й?, I = 1,2,3, р" = ^ р", и" = \ г| I (4)
для узловых перемещений и и давления р получим результирующую связанную систему уравнений на каждом конечном элементе Vei в каждый момент времени tm:
fL Dl-(®lm = -( 0 (5)
ID tK! Ipj {Du + nßp + tF) k >
где:
wfc
^дгайик : С : дгай\\?к dУ = J Сцыи^Ющ йУ — Ь^й?
Те1 Уе1
Уе1 уе1
^дгайрк ■ итм йУ = | р^м^ = £>£чррй?
1Шт»Ч» сШ = и?лЧ» сШ = (6)
ПР ^ Г „ ПР
Г ПР Г N N ПР I (дгайр • к • дгайц--рц) йУ — I (к^щц"--рц) йУ — — Крчррцч,
Те1 Уе1
крч = | (кцМрЛМчЛ - ^ ИрИч)йУ
I вцйЕ = I йЕ = рр= I СМР° йЕ
Как было отмечено выше, для реальных многослойных грунтов матрица жесткости системы, получаемая ассемблированием (6), является плохо обусловленной. Это, по-видимому, связано также с тем, что при малом т матрица системы (5) «похожа» на матрицу с нулевым углом, типичную для задач несжимаемой упругости. Как известно [9], в последнем случае аппроксимация давления р должна осуществляться способом, отличным от три-линейных базисных функции, использованных в (6). Вычислительные эксперименты показывают, что обращение матрицы (5) неустойчиво. Для разрешения проблемы можно использовать другой вариант аппроксимации по времени.
3. Во-первых, можно использовать неявную схему, в которой оба уравнения имеют первый порядок аппроксимации 0(т) относительно шага т. Для удобства читателя, дальнейшее изложение будем вести, используя обозначения дифференциальных операторов вместо подматриц (5). Отнесем левую часть первого уравнения (1) к временному слою tm, а правую к £т_1. Второе уравнение в (1) запишем на временном слое £т_х. Тогда конечномерная аппроксимация уравнений (1), а точнее говоря уравнений (3), примет вид:
—сИу(С\ дгасС) О
О —тсИу(к ■ дгай) + пр
( ит ]
• [рт-1 ]
( дгайрт~г |
Шуи™-1 - й(уйт~2 - пррт"2 - тFm-1J ( )
Такое разделение уравнений по времени позволяет решать их независимо одно от другого. Теперь устойчивость решения линейной системы (7) на каждом временном шаге определяется для первого уравнения только разбросом компонент тензора С, а для вто-
6/2П11 ВЕСТНИК
_6/2011_МГСУ
poro, соответственно тензора fe. Отношение этих параметров между собой и наличие множителя т в одном из блоков матрицы на устойчивость обращения не влияют.
Во-вторых, можно использовать неявную схему с внутренними итерациями, предложенную еще в [2]. Для решения системы (5) применим итерационный процесс с внутренними итерациями на каждом шаге по времени:
( div(C: gradum-s) — gradpm-s~1 — О
l-divüm-s~1 + Tdiv(k ■ gradpm-s) - n^pm-s~1 - ( )
= -(divum~1 + nppm~1 + rFm~1)
где s — 1,2,3,... - номер внутренней итерации. Данный итерационный метод (8) также позволяет «развязать» систему уравнений Био. Как показано в работе [2], скорость сходимости такого итерационного метода не зависит от шагов сетки, а следовательно, числа уравнений решаемой линейной системы. Также важно, что она не зависит от разброса значений самого «дисперсного» коэффициента задачи - коэффициента фильтрации грунта fe, но зависит от не столь контрастных коэффициентов - тензора модулей упругости грунта и пористости. Это увеличивает шансы на устойчивое и надежное решение системы (5).
4. Авторами написан пакет программ на языке FORTRAN, реализующий приведенные выше алгоритмы в трехмерной постановке. Ниже приведем в качестве примера численный расчет деформаций грунта и изменения порового давления нефти из месторождения, расположенного в Западной Сибири. Месторождение моделируется шес-тислойной средой с горизонтальным расположением слоев. Данные, с которыми производились вычисления, приведены в таблице 1. Откачка осуществляется из двух скважин, находящихся в верхних 50 метрах в 5-м слое, суммарным расходом G=3000 м3/сут. Расстояние между ними 2 км. Протяженность моделируемого месторождения следующая: длина 10 км, ширина 8 км, глубина 8 км. Граничные условия моделируют непроницаемость месторождения по внешней границе. Результаты расчета приведены в рисунках 1 - 4.
Таблица 1.
№ слоя Мощность, м Модуль упругости, кг/см2 Коэффициент Пуассона Коэффициент фильтрации, м/сут Пористость, %
1 200 5000 0,25 5 25
2 400 90000 0,35 0,000001 40
3 800 130000 0,28 0,1 20
4 900 120000 0,3 0,000001 15
5 200 160000 0,28 0,1 20
6 5500 200000 0,25 0,000001 5
Z Д
\
% V \ E -0.04362 £ l 1
\ ™ 0> " ""¡„„„ena».»^'"
Рис.1 Осадка земной поверхности Рис. 2 Поровое давление в среднем
месторождения. сечении в 3-м слое.
Рис. 3 Горизонтальные перемещения иу Рис. 4 Напряжения сдвига ayz скеле-скелета в среднем сечении в 5-м слое. та в среднем сечении в 4-м слое.
Работа выполнена при поддержке грантов РФФИ № 08-05-00578-а, 07-01-92111-
ГФЕН_а.
Литература
1. Терцаги К. Теория механики грунтов. - М.: Госстройиздат, 1961. - 507 с.
2. Шешенин С.В., Киселев Ф.Б., Численное моделирование нестационарной фильтрации в грунте //Изв. РАН, МТТ, 3, 1996.
3. Biot M.A. Theory of elasticity and consolidation for a porous anisotropic solid. // Journal of Applied Physics. 1955. Vol. 26. No. 2. P. 182-185.
4. Gambolati G. Numerical models in land subsidence control. // Computer Methods in Applied Mechanics and Engineering. 1975. No. 5. P. 227-237.
5. Jacob C.E. Flow in groundwater. // Rouse h. engineering hydraulics. New York: Wiley, 1950. P. 321-386.
6. Kalinin E.V., Sheshenin S.V., Artamonova N.B., Kiselev F.B. Numerical investigations of the influence of fluid extraction upon the stress state of the rock masses. // Eng. Geology and Environment. Mat. Intern. Symp. Athens, Greece, 1997. P. 725-728.
6/2Q11 мвВЕСТНИК
7. Tran D., Settari A., Ngheim L. New Iterative Coupling between a Reservoir Simulation and a Geomechanics Module. SPE Journal, Vol. 9, No 3, Sept. 2004, P. 362-369.
8. Wan J., Durlofsky L.J., Hughes T.J.R., Aziz K. Stabilized Finite Element Methods for Coupled Geomechanics - Reservoir Flow Simulations, SPE Paper 79694, SPE Reservoir Simulation Symposium, 3-5 February 2003, Houston, Texas.
9. Zienkiewicz O.C. and Taylor R.L., The Finite Element Method, 5th edition Vol. 2: Solid Mechanics, Butterworth-Heinemann, 2000.
References
1. Biot M.A. Theory of elasticity and consolidation for a porous anisotropic solid. // Journal of Applied Physics. - 1955. - Vol. 26. - No. 2. - P. 182-185.
2. Tercagi K. Teorija mekhaniki gruntov. - M.: Gosstrojjizdat, 1961. - 507 s.
3. Jacob C.E. Flow in groundwater. // Rouse h. engineering hydraulics. - New York: Wiley, 1950. - P. 321-386.
4. Gambolati G. Numerical models in land subsidence control. // Computer Methods in Applied Mechanics and Engineering. - 1975. No. 5. - P. 227-237.
5. Sheshenin S.V., Kiseljov F.B., Chislennoe modelirovanie nestacionarnojj fil'tracii v grunte //Izv. RAN, MTT, 3, 1996.
6. Kalinin E.V., Sheshenin S.V., Artamonova N.B., Kiselev F.B. Numerical investigations of the influence of fluid extraction upon the stress state of the rock masses. // Eng. Geology and Environment. Mat. Intern. Symp. Athens, Greece, 1997. - P. 725-728.
7. Wan J., Durlofsky L.J., Hughes T.J.R., Aziz K. Stabilized Finite Element Methods for Coupled Geomechanics - Reservoir Flow Simulations, SPE Paper 79694, SPE Reservoir Simulation Symposium, 3-5 February 2003, Houston, Texas.
8. Tran D., Settari A., Ngheim L. New Iterative Coupling between a Reservoir Simulation and a Geomechanics Module. SPE Journal, Vol. 9, No 3, Sept. 2004, pp. 362-369.
9. Zienkiewicz O.C. and Taylor R.L., The Finite Element Method, 5th edition Vol. 2: Solid Mechanics, Butterworth-Heinemann, 2000.
Ключевые слова: модель Био, трехмерная постановка, связанная задача, слоистая среда, плохая обусловленность, неявная схема, внутренние итерации.
Key words: Biot's model, three-dimensional problem, coupled equations, layered medium, ill-conditioning, implicit difference scheme, internal iterations.
Авторы:
1. Шешенин Сергей Владимирович, д.ф.-м.н., профессор (МГУ им. М.В.Ломоносова), 117234, РФ, г. Москва, Ленинские горы, д. 1, тел./факс: +7(495) 939-43-43, e-mail: sheshe-
ni@mech.math.msu.su
2. Киселев Федор Борисович, к.ф.-м.н., доцент (ГОУ ВПО МГСУ), 129337, РФ, г. Москва, Ярославское шоссе, д. 26, тел./факс: +7(499) 183-24-01, e-mail: afi-ko@mail.ru 3. Артамонова Нина Брониславовна, канд.геол.-мин.наук, с.н.с. (МГУ им. М.В.Ломоносова), 117234, РФ, г. Москва, Ленинские горы, д. 1, тел./факс: +7(495) 939-43-43, e-mail:
artamonovanb@mail.ru
Рецензент: Звягин Александр Васильевич, д.ф.-м.н., профессор кафедры волновой и газовой динамики МГУ им. М.В.Ломоносова
