CC BY
54
УДК 539.3
О ЧИСЛЕННОЙ РЕАЛИЗАЦИИ МОДЕЛИ БИО.
Шешенин С.В., Киселев Ф.Б. (Москва)
Abstract
Some measurements of pressure in observation wells situated near a pumping show well the increase of pressure instead of predicted decrease. The approach based on simplified model of consolidation fails to describe anomalous pressure behavior. More complicated Biot's model has been taken into our consideration.
Finite element method is used for spatiai digitization. Time discretization and internal iterations at each time step is used.
Одной из задач механики гетерогенных сред является задача о движении жидкости в порах грунта под действием массовых сил и гидродинамического давления. С фильтрацией имеют дело при описании притока жидкости к дренажным щелям и буровым скважинам, просачивания воды под плотинами, таяния ледников, снежных пластов и т.д. В предлагаемой статье используется модель Био [1] для численного моделирования нестационарной фильтрации сжимаемой жидкости в грунтах и горных породах, моделируемых упругой пористой средой, обычно именуемой скелетом (каркасом). Стационарная (установившаяся) фильтрация хорошо описывается упрощенной моделью, называемой в литературе моделью упругого режима. Эта модель вытекает при некоторых предположениях из модели Био и описывается несвязанными уравнениями упругого равновесия каркаса и диффузии жидкости. В случае нестационарной фильтрации, например, в начале откачки из скважины или в случае поднятия воды в шлюзах эта модель может приводить даже к качественным расхождениям с данными измерений [2,3]. В таких случаях для расчета фильтрации следует применять полную модель Био, состоящую из связанных уравнений равновесия упругого каркаса и диффузии жидкости. Однако численная реализация этой модели приводит к трудностям, один вариант преодоления которых описывается в данной статье.
1. Первое уравнение модели Био имеет вид
Lu =div C: grad u =grad p, (1)
где u - вектор перемещения каркаса, C - эффективный тензор модулей упругости пористого грунта, p - давление жидкости. Здесь предположено, что грунт полностью насыщен жидкостью, а давление, перемещения скелета и другие далее вводимые величины представляют собой разности этих величин в состоянии движения и состоянии покоя под действием массовых сил. Поэтому последние исключены из рассмотрения.
Второе дифференциальное уравнение, вытекающее из уранения неразрывности и закона Дарси, записывается в виде
div ( • grad p) = divU + npp . (2)
При выводе данного уравнения предполагается, что жидкость баротропна и в - ее сжимаемость, а материал скелета несжимаем. Здесь к - тензор проницаемости пористой среды, П - ее пористость. В этом уравнении также произведена замена полной производной от давления по времени на частную, т.к. предполагается малость скорости фильтрации . Тензоры С и к мы предполагаем ортотропными.
Таким образом, система уравнений (1),(2) описывает фильтрацию жидкости и вызываемую последней деформацию скелета. Будем рассматривать краевую задачу для системы (1),(2) со следуюшими краевыми условиями, соответствующими задаче об откачке жидкости через скважину (рисунок):
на
2^ : V'С : gradu=0, р =0,
(3)
на
^2 : и = 0,
дР
дv
=0,
(4)
на :и -у = 0, т-а-у = 0, г^г —— = 1/2п1 Q,
3 ду
на
^4 : и - у=0, т -а - у=0, др =0.
ду
(5)
(6)
Рис. Схема к постановке задачи
В последних формулах через I обозначена длина скважины, V, Т - соответственно нормальный и касательный единичные вектора к границе, Q - объем откачиваемой (закачиваемой) жидкости в единицу времени. Начальные условия имеют простой вид:
u=0, p=0 при t = 0. (7)
Итак, краевая задача описывается уравнениями (1),(2) с краевыми условиями (3)-(6) и начальными условиями (7).
2. Для исследования описанной выше краевой задачи сформулируем
соответствующую ей вариационную задачу. Для этого умножим уравнение (1) на функцию W, а уравнение (2) на функцию q, проинтегрируем по объему области и, учитывая граничные условия (3)-(6), получим:
J grad u : C: grad w dV + J"^ grad p w dV =0, (8)
grad p ■ k • grad q dV +(diva + nfip) qdV = 1/2П Je QdE . (9)
Система вариационных уравнений (8),(9) определяет обобщенное решение исходной краевой задачи как пару функций и,р, принадлежащих соответствующим
подпространствам Ни и Н пространства функций Ж2' , определяемых ограничениями,
накладываемыми граничными условиями (3)-(6).
В результате конечноэлементной аппроксимации по пространственным координатам и замены производной по времени разностной производной, вместо (8) можно получить систему вида
Г A i 2 4 S u Г-J 0
I 21Л тЛ22 + nP_ Ip! =1 Л um2 + nppm2 + TfQ W
T T
(причем A12 = A21 ), где Т - шаг дискретизации по времени, а mx = m -1,m2 = m для явной схемы и mx = m, m2 = m — 1 для неявной. Использование явной схемы невозможно по ряду причин. Во-первых, ее можно реализовать только для несжимаемой среды, когда в= 0 . Но в этом случае она неустойчива. Кроме этого, такая схема только условно может быть названа явной, поскольку требует обращения несимметричной матрицы A21. Использование же неявной схемы приводит к необходимости обращения блочной матрицы, стоящей в квадратных скобках. При малом шаге по времени эта матрица плохо обусловлена. Кроме того, вряд ли можно ожидать ее знакоопределенности. Ниже мы опишем способ преодоления указанных математических трудностей.
Решение вариационного уравнения (8) при фиксированной функции grad p существует и единственно [4,5]. Следовательно, определен оператор L \ действующий из ь2 в W1, такой, что U = L 1 ■ grad p. Обозначение вызвано тем,что через L мы обозначили оператор теории упругости (1), к которому данный оператор является обратным. Подставляя последнее выражение в (9),получим вариационное уравнение относительно функции p:
j gradp ■ к ■ gradqdV +J [divL1 gradp + nfip)dV =
= 1/2nl Je QdE. (10)
От уравнения (10) можно стандартным образом перейти к операторному уравнению в некотором гильбертовом пространстве HD [4,5], которым может служить
подпространство W2, (учитывающее граничные условия на функции u, p), или "энергетическое" пространство HD, состоящее из того же множества функций, но с другим скалярным произведением. В итоге получаем операторное эволюционное уравнение
Bp = Ap + f, (11)
где операторы A, B и "правая" часть f определяются соотношениями:
(Ap,q)D = — JVgradp■ k ■ gradqdV ,
(B p, q)D = JV (divL1 gradp+nfi p) qdV , (12)
(f,q)D =12 nl JE3QqdV .
Из (12) видно, что A является отрицательно определенным оператором. Можно показать, что в силу граничных условий (3)-(6) оператор B также является симметричным и положительным. Поэтому определен B 1 [6] и уравнение (11) может быть переписано в виде p = B-Ap + B 1 Af, откуда сразу следует, что оно эволюционного типа и, следовательно, исходная краевая задача корректна. Рассмотренное преобразование исходной краевой задачи к уравнению (11), естественно, можно осуществить и в конечномерном случае, произведя пространственную дискретизацию, например, методом конечных элементов.
3. В данном пункте мы рассмотрим неявную разностную схему для решения уравнения (11), записываемую в виде
pm1 — pm—1
B£------I—=Apm + f, (13)
Т
где tm = Tm, m = 0,1,2... Хорошо известно [ 6 ], что схема (13) абсолютно устойчива в пространстве H A как по начальным данным, так и по правой части. Явную схему использовать нецелесообразно из-за слишком мелкого шага по времени и трудностей при обращении оператора B.
Обращение на каждом временном шаге оператора B — tA нисколько не сложнее обращения оператора B , хотя и является нетривиальной задачей, для решения которой мы используем внутренний итерационный процесс. Для того чтобы описать данный процесс, перепишем (13) в виде
Aopn + F = 0, (14)
A0 = A — 1 tB, F = 1/ tB pm—1 + f. y ’
Для решения (14) применим стандартный двухслойный итерационный процесс
рШ^+1 _ рШ.в
В^---------+ А0 рт*+^=0, (15)
а
сходимость которого будет иметь место при а<2 у2 , если выполняются неравенства энергетической эквивалентности операторов А0 и В0 [ 4, 5, 6 ]:
0 </,( _ В„)<(_А,)<у2 (_ Во). (16)
Можно получить оценки для констант, входящих в неравенства (16). Прежде всего выберем оператор В0 = А _ аЕ, где Е - единичный оператор, а - параметр, выбор которого описан ниже. Для удобства проведения оценок введем операторы С и С0 :
(С р, д)в = - V &таёд ■ А1 ^уаёрёУ ,
(Со р, д)0 = — У $гаёд ■ ( Ь1 ■ ^аёр) ёУ .
(17)
Поскольку оператор Ь спектрально эквивалентен оператору Лапласа А с компонентами А. = А 8. , то имеют место неравенства 0 < С1 (—А) < (—Ь) < С2 (—А) . Учитывая (17), будем иметь
1/с Со <С Со. (18)
Относительно собственных значений оператора С0 известно, что они расположены на отрезке [ т, 1], где т > 0 . Поскольку В = С + П^Е и т/С2 <С < 1/С2 Е , то
(тС + пв) (—Е) <— В < (1/с + пв) (—Е). (19)
Выбирая а=1/ Т, получим :
У\ = (ЧС2 + пв) , 7 2 = (1/С1 + пв) . (20)
Скорость сходимости процесса (15) определяется отношением ух1 у2 и, следовательно, не
зависит от проницаемости к , но зависит от значений упругих модулей С. Если коэффициенты С.к1 удовлетворяют обычным неравенствам / к.к. <С.к1 И/ < к2 ИМ. , то
С = МИ 1, С2 = /^ , где М - константа Корна. Поскольку в у1, у2 имеются две константы
т и М , значения которых неизвестны, то целесообразно применять адаптацию итерационного параметра [ 7 ] в ходе итераций (15) .
На практике итерационная схема (15) применяется к решению конечномерной задачи. Достоинством предложенного метода является то, что вследствие спектральной эквивалентности разностных аналогов операторов А0 и В0 обеспечивается сходимость итераций со скоростью, не зависящей от числа уравнений. Сама же неявная схема безусловно устойчива, что оставляет свободу в выборе шага по времени.
Литература.
1. Biot M.A. General theory of three-dimensional consolidation// Journal of applied physics.-1941, Vol. 12, P. 155-164.
2. Wolff R.G. Relationship between horisontal strain near a well and reverse water level
fluctuation // Water resources research.- 1970, Vol. 6, No 6, P. 1721-1728.
3. Borja R.I. One-Step and Linear Multistep Methods for Nonlinear Сonsolidation //
Computer Methods in Applied Mechanics and Engineering.- 1991, Vol. 85, No 3, P. 239-
272
4. Дьяконов Е.Г. Минимизация вычислительной работы. Асимптотически оптимальные алгоритмы для эллиптических задач.- М. : Наука, 1989.- 272 с.
5. Победря Б.Е., Шешенин С.В., Холматов Т. Задача в напряжениях.- Ташкент: ФАН, 1988.- 200 с.
6. Самарский А.А. Теория разностных схем.- М.: Наука, 1989.- 616 с.
7. Хейгеман Л., Янг Д. Прикладные итерационные методы.- М.: Мир, 1986.-448 с.
