Содержание к диссертации
Введение
1 Постановка задач теории упругости. несвязанные модели био процессов фильтрации в пористых и трещиновато-пористых средах 13
1.1 Математическая модель распространения волн в упругой среде 13
1.2 Модель Био 15
1.3 Обобщение модели Био на случай трещиновато-пористой среды 18
2 Конечно-разностная модель на неортогональных сетках для эллиптического уравнения с разрывным коэффициентом. метод опорных операторов. сходимость . 20
2.1 Интегральные тождества с методе опорных операторов с разрывными коэффициентами 20
2.2 Определение метрического тензора 23
2.3 Аппроксимация потоков и сильная сходимость 26
2.4 Метод опорных операторов в задаче теории упругости с разрывными коэффициентами 29
3 Конечно-разностная модель на неортого нальных сетках для параболического уравнения и для уравнений фильтрации в задачах с двойной пористостью. устойчивость.сходимость . 33
3.1 Слабосходящийся алгоритм решения параболического урав-ненияс разрывными коэффициентами . 33
3.2 Регуляризованная схема 39
3.3 Метод опорных операторов в задаче теории фильтарции . 45
4 Алгоритмы осреднения в ячейках для эл липтического скалярного уравнения и системы уравнений теории упругости. слабая сходимость . 51
4.1 Об алгоритмах осреднения в ячейках для эллиптических кра евых задачах 51
4.2 Обобщение схемы averaging'a 54
4.3 Метод averaging'a в задаче теории упругости 56
4.4 Трехмерная полностью анизотропная задача теории упругости 61
5 Результаты расчетов 66
Заключение 79
Литература 81
- Обобщение модели Био на случай трещиновато-пористой среды
- Метод опорных операторов в задаче теории упругости с разрывными коэффициентами
- Метод опорных операторов в задаче теории фильтарции
- Трехмерная полностью анизотропная задача теории упругости
Введение к работе
Математические модели в подземной гидродинамике, описывающие влияние потока флюида на напряженно-деформированное состояние окружающей среды, всегда представляли практический интерес. Например, при разработке гигантских углеводородных месторождений, характерных для России, возможно влияние процессов фильтрации на сейсмическую активность региона. При этом потоки флюида приводят к возникновению отсутствующих в геостатике касательных напряжений [Турунтаев и др., 1994].
Впервые подобная модель была предложена Терцаги [Terzaghi, 1936]. Он ввел эффективный тензор напряжений, зависящий от деформации скелета и смещений флюида. Однако применение соотношений Терцаги на практике представляет большие сложности. Реально могут быть наблюдаемы деформации порового объема как целого, но не деформации скелета и смещения флюида по отдельности. Сейчас процессы фильтрации в насыщенных пористых средах, сопровождаемые изменением напряженно-деформированного состояния, обычно рассматриваются в рамках квазистатической модели Био [Biot, 1962; Rice and Cleary, 1976]. В зависимости от того, рассматривается ли обратное влияние изменения напряженно-деформированного состояния на движение флюида, модель называется связанной или несвязанной.
В работах [Авербух и др., 1995; Гнедин и др., 1996] в рамках несвязанной модели Био были получены пространственные распределения давления фильтрующегося флюида, которые затем использовались для определения полей деформации и напряжений. В отличие от [Desrocher, 1994], где определялась только величина относительного изменения объема, в указанных выше работах были вычислены все компоненты тензора напряжений. В работе [А.В. Колдоба и др., 1999] была рассмотрена двумерная связанная задача Био, т.е. совместное решение уравнение фильтрации и уравнений теории упругости.
Однако в моделях такого типа пористая среда рассматривается как совокупность непроницаемых зерен, разделенных порами. В действительности же для всех естественных пластов характерна развитая в той или иной степени трещиноватость. Здесь возможны два подхода. В первом из них рассматривается более или менее регулярная система трещин. Однако, даже предположив, что мы сможем решить задачу с достаточно общим их распределением, все равно, находясь на поверхности Земли, невозможно
достоверно судить о конфигурации этой системы. Здесь ситуация аналогична введению модели пористой среды. Даже если бы мы смогли проинтегрировать уравнения движения вязкой жидкости в порах, все равно их расположение в среде нам неизвестно. Так были введены осредненные величины, такие как пористость, давление, проницаемость, скорость фильтрации, а все основные законы формулировались в терминах этих величин. Именно по этому пути пошли сторонники второго подхода. Так была введена модель фильтрации в трещиноватых средах [Баренблатт и др., 1960].
В связи с вышесказанным представляет интерес построение единой модели фильтрации в трещиноватых средах и эволюции напряженно-деформированного состояния среды.
При этом увеличение давления приводит к появлению сдвиговых деформаций, которые могут стать причиной оползания пород на границе разлома. Такой процесс может спровоцировать сейсмическую активность среды вплоть до техногенного землетрясения.
Для описания фильтрации в трещиноватых средах используется модель с двойной пористостью и двойной проницаемостью [Баренблатт и др-, 1960]. Насыщенная среда представляет собой структуру, составленную из пористых блоков и трещин. Если макроскопическое описание пористой среды есть результат усреднения по масштабам большим размеров пор, то для трещиноватых сред вводятся характеристики, являющиеся результатом усреднения по масштабам большим сравнительно с размерами блоков. Отличие развиваемой здесь схемы от обычной схемы фильтрации в пористой среде состоит во введении в каждой точке пространства двух давлений жидкости - давления жидкости в порах и давления жидкости в трещинах -и учета обмена жидкостью между трещинами и порами. При определенных предположениях получается выражение для интенсивности этого обмена. Выводятся основные уравнения фильтрации жидкости в трещиноватой породе и более общие уравнения фильтрации жидкости в пористой среде с двойной пористостью.
Таким образом, уравнения фильтрации в трещиновато-пористой среде представляют собой уравнения с двойной пористостью и двойной проницаемостью. Причем, поскольку проницаемость в трещинах очень велика, давление в трещинах устанавливается гораздо быстрее, чем в порах, так что для его описания может быть использовано квазистатическое уравнение. Кроме того, можно предположить, что объем трещин пренебрежительно мал, и, как следствие, можно пренебречь величиной флюидосодержания в трещинах. В результате получается модель, отличающаяся от модели Био
для пористых сред уравнением для давления в трещинах при условии, что обобщение закона Гука в модели Био сохраняет свое выражение. В дальнейшем мы используем несвязанный вариант данной модели, предложенный Николаевским [1996].
Хорошо известно, что для земной коры характерна структурная рас-слоенность. Эта расслоенность проявляется в чередовании зон высоких и низких скоростей, областей повышенных и пониженных напряжений. Существует множество примеров предельных случаев, когда фильтрация происходит в достаточно узких областях - разломах.
Необходимость учета геологической структуры среды, которая состоит из ряда слоев сложной геометрии, резко отличающихся по своим свойствам, привела к развитию численных методов решения разностных задач на криволинейных сетках, адаптированных к структуре среды. Построение конечномерной аппроксимации задачи на криволинейных сетках может быть осуществлено как проекционными методами, например, методом конечных элементов, так и методом конечных разностей. Если использовать метод конечных элементов с треугольными элементами, то в двумерном случае могут быть получены пятиточечные схемы для ортогональных сеток или семиточечные для неортогональных сеток при соответствующей ориентации прямоугольников. При соответствующем способе аппроксимации источников в такой задаче может быть доказана сходимость алгоритма в метрике Н1. В настоящей работе рассмотрен для данной задачи метод опорных операторов. Особенностью этого метода является использование при определении разностных операторов дивергенции и градиента с помощью разностного аналога известного интегрального тождества
/ р- VudV + / udivpdV — I pnudS (*)
Ja Jn Jm
To есть один этих из операторов определяется вручную (он называется опорным), а другой находится как сопряженный к нему.
Можно показать, что данный метод консервативен. Кроме того, разностная схема, построенная этим методом на четырехугольных сетках, переходит в обычную пятиточечную, если все ячейки прямоугольные, но, в отличие от метода конечных элементов на треугольных ячейках, имеет девятиточечный шаблон в общем случае.
Как для скалярных эллиптических уравнений, так и для стационарных задач теории упругости в двумерном случае метод опорных операторов достаточно подробно исследован, однако только для постоянных коэффици-
ентов [Самарский и др., 1996]. В этой работе доказана сходимость первого порядка как в энергетической, так и в р-нормах.
Основная проблема для построения разностных схем для уравнения теплопроводности - это вычисление нормального потока через грани ячейки. Эта проблема распадается на две задачи: определение компонентов потока и вычисление баланса тепла, что позволяет определить выражение для дивергенции вектора теплового потока. При этом для определения величины потока через каждую грань могут быть использованы как ковариантные, так и контравариантные компоненты вектора потока. Первые определяются очевидным образом как разностные производные вдоль ребер сетки. При вычислении потока через грань необходимо уметь вычислять контравариантные компоненты вектора потока. В работе применяется предложенный А.А. Самарским и др. [1996] метод, основанный на построении специального метрического тензора. После этого возможно вычисление балансов для каждой ячейки и определение сопряженного оператора в силу (*) Особенностью предложенного метода является девятиточечные схемы, которые обладают более высокой точностью по сравнению с семиточечными схемами метода конечных элементов.
На основе метода опорных операторов были построены экономичные абсолютно устойчивые алгоритмы решения параболических уравнений с разрывными коэффициентами [Zaslavsky, 2003]. Экономичность схемы означает, что число операций для получения решения на неявном слое по времени пропорционально числу узлов сетки. Очевидно, что явная схема экономична. Однако существует жесткое ограничение на шаг по времени для устойчивости схемы: т ~ h2. При этом полностью неявная схема абсолютно устойчива, но не экономична. Кроме того, она является схемой первого порядка аппроксимации по времени. В соответствии с указанными требованиями были построены регуляризованная схема, а также схемы на основе разностного и вариационно-разностного метода потоков. Кроме того, были найдены условия абсолютной устойчивости для регуляризованной схемы.
Сейчас метод опорных операторов широко используется, в том числе и в прикладных задачах. В частности, на его основе построен разностный алгоритм решения связанной задачи пороупругости в плоской геометрии [Колдоба и др., 1999].
При решении квазистатической задачи Био, а именно ее упругой части, избежать обращения матрицы не удается. Хорошо известно, что при использовании любого итерационного алгоритма решения линейной системы точность полученного решения после заданного числа итераций зависит от
числа обусловленности обращаемой матрицы. Иначе говоря, чем оно больше, тем больше число итераций необходимо сделать, чтобы получить решение с данной точностью. Но если разлом достаточно узок по сравнению с характерным размером всей области, адаптированная сетка будет достаточно подробной, что ведет к большим затратам процессорного времени для решения системы.
В связи с этим возникает потребность считать на более крупных ячейках - с размерами, большими, чем ширина разлома. Но в этом случае некоторые ячейки будут достаточно сложной структуры, состоять из разных материалов. Таким образом, необходимо уметь вычислять осредненные коэффициенты в этой ячейке таким образом, чтобы полученная математическая модель учитывала свойства начальной среды.
Проблема осреднения, или upscaling'a или averaging^, возникла еще при осреднении коэффициента проницаемости в задачах фильтрации. Хорошо известно, что в одномерном случае эффективный коэффициент в ячейке от Х{ до Xi+\ есть среднее гармоническое
h h
Jii к
где h = хі+і —Xi. Такое среднее значение дает первый порядок сходимости в энергетической норме полученной разностной схемы.
Перейдем к рассмотрению двумерной задачи. Принципиальной особенностью метода осреднения в этом случае является изменение характера среды. Действительно, если изначально среда в каждой ячейке изотропна, то малый, но конечный фрагмент среды, содержащий структурные особенности, не обладает свойством изотропии. Например, если линии разрыва параллельны одной из осей, например у, то в направлении оси х используется среднее гармоническое, а по оси у - среднее арифметическое
jHkds
где Н - ячейка, по которой производится осреднение. В общем случае расположения разрывов наиболее часто при решении таких задач используется, в том числе и в прикладных пакетах, осреднение Кинга [King et al., 1993], которое было построено специально для изотропной среды, причем двумерная область состоит из прямоугольников, отличающихся по структуре. Автор, проводя аналогию между законом Дарси и законом Ома, находит эффективное сопротивление в электрической сети. Показано, что
алгоритм Кинга равносилен определению эффективного коэффициента из соответствующей простейшей разностной схемы [Максимов и др., 2002]. Однако эффективный тензор, определенный по указанному алгоритму, диагоналей. То есть, в двумерном случае при осреднении коэффициента теплопроводности к возможно определить лишь две константы кх и ку. При этом очевидно, что градиент температуры вдоль одной из осей создает в любой малом, но конечном элементе среды сложной структуры тепловой поток по обеим осям в общем случае. Таким образом, в эффективном тензоре должны присутствовать и недиагональные компоненты.
Эффективный тензор такого вида получен на Лебедевской сетке Друс-киным [Moskow, Druskin et al., 1999]. Он рассматривал эффективный тензор как матрицу квадратичной формы, соответствующей энергии системы. В стандартном averaging'e эффективный тензор в базисе (т; п), связанном с разрывом (ft нормален к нему), имеет вид
(ІгО-
В узловом averaging'e он определяет эффективный тензор, исходя из принципа совпадения разностной и непрерывной энергий на любом элементе линейной оболочки, натянутой на некоторые базисные вектора, при этом решение может быть аппроксимировано элементом этой линейной оболочки. Из физических соображений очевидно, что принцип тождества энергий - один из самых важных для совпадения некоторых свойств исходной среды и модели после осреднения. Аналитически же удалось доказать слабую (в смысле энергетического скалярного произведения) сходимость с первым порядком решения полученной разностной задачи к точному решению исходной задачи.
Таким образом, в результате осреднения элемент изначально изотропной среды, но имеющий сложную внутреннюю структуру, в общем случае может быть адекватно описан только с помощью полностью анизотропной модели.
Особый интерес представляет построение алгоритмов осреднения в ячейках для динамических задач теории упругости. Во избежание использования сеток, адаптированных к структуре среды, построение которых в общем случае достаточно сложно, возможно проводить расчеты на прямоугольных сетках, предварительно проводя осреднение коэффициентов в них. Как один из способов решения обратных задач, то есть определения параметров среды, часто используется метод подбора. Поэтому необходимо
точно определять эффективные тензоры в ячейках. От этого зависит точность воспроизведения упругих волн, в частности, волн Стоунли. Последние как раз играют важнейшую роль в задачах, связанных с акустическим каротажем.
Задачи теории упругости с точки зрения принципа тождества энергий были рассмотрены в плоской геометрии [В.П. Мясников, М.Ю. Заславский, А.Х. Пергамент, 2004]. Легко видеть, что и в теории упругости эффективный тензор является в общем случае также полностью анизотропным.
Задачи пороупругости с этой точки зрения являются более сложными в том смысле, что осреднять необходимо не только коэффициенты Ламе (в изотропном случае), но и источники, обусловленные потоком флюида.
Таким образом, рассматриваемые проблемы включают решение параболических уравнений в областях сложной формы и решение эллиптических краевых задач с источниками, характерные масштабы которых много меньше актуальных размеров задачи. Для решения таких задач обычно используются сетки, адаптированные к структуре среды, а из численных методов - методы конечных элементов, конечных объемов (I.Babuska et al., 1994; Glowinski R., Wheeler M.F., 1988), а в отечественных исследованиях - метод опорных операторов. Следует отметить, что если разностные операторы градиента и дивергенции являются сопряженными друг другу, а также имеется аппроксимация интеграла энергии, то разностная схема является консервативной, и можно установить слабую сходимость алгоритма. Кроме того, если имеет место аппроксимация потоков, то можно доказать сильную сходимость в энергетической норме и, как следствие, в силу теорем вложения, в норме Lp.
Как правило, задачи аппроксимации решаются для линейных подпространств и многообразий, причем базисные функции этих множеств должны быть таковы, чтобы решение задачи могло бы быть аппроксимировано элементами линейной оболочки с достаточной точностью.
Далее требуется равенства непрерывной и дискретной энергий на элементах линейного подпространства или равенства непрерывного и разностного потоков через грани ячеек для тех же функций. Первый факт ведет с слабой сходимости, из второго же требования можно вывести сильную сходимость алгоритма.
Из вышеприведенных рассуждений, в частности, следует, что метод конечных объемов консервативен, но обладает слабой сходимостью, так как для четырехугольной ячейки (или шестигранной в трехмерном случае) обычно не обеспечивается равенство непрерывного и разностного потоков
для элементов span{l, х, у} (или span{l, зг, у, z} в трехмерном случае). Что касается метода конечных элементов, то при использовании ячеек в виде треугольников (или тетраэдров) удается построить алгоритм, который обладает сильной сходимостью. Увеличение числа узлов или граней требует расширения базиса, точнее использования полилинейных аппроксимаций на основе элементов span{l,x,y,xy].
Метод опорных операторов в отличие от метода конечных объемов позволяет построить алгоритм, сходящийся в энергетической норме и в норме Ьр для двумерных задач. Авторами метода построены аппроксимации для скалярных произведений, порожденных интегралами энергии, на основе введения своеобразного метрического тензора, определенного в ячейках двумерной сетки. При этом, в согласии с вышеизложенной схемой рассуждений выражение для скалярного произведения (и, соответственно, интеграла энергии) является точным для элементов линейной оболочки, натянутой на совокупность базисных векторов. Выражение для метрического тензора содержит неизвестные параметры, выбор которых позволяет обеспечить аппроксимацию потоков.
Особенностью метода опорных операторов является построение сопряженного объема, для которого выполняются условия, обеспечивающие аппроксимацию потоков. Следует отметить, что в трехмерном случае при наличии шестигранных ячеек не удается обеспечить и аппроксимацию энергий и явным образом построить сопряженные объемы. Однако на сетке из тетраэдров эти объемы можно явно найти, и искомый алгоритм будет схож с построенным в двумерном случае.
Авторы многочисленных работ как по методу конечных объемов, так и по методу опорных операторов фактически использовали в качестве базисов линейные и полилинейные функции. Это корректно, если решение не имеет особенностей. Если же коэффициенты уравнения разрывны, то необходимо построить базисы, позволяющие учитывать особенности решения. В качестве таких базисов были выбраны предложенные I.Babuska функции, которые в двумерном случае сводятся к span{l,m ^,/J1'7"^}. Здесь г - радиус-вектор, {т, п} - касательный и нормальный к разрыву вектора, образующие правый базис, к - коэффициент уравнения, в окрестности разрыва зависящий только от ft г. В результате, используя вышеизложенную схему, удалось построить для двумерных задач фильтрации в пористо-трещиноватых средах и задач теории упругости упругости алгоритм метода опорных операторов, сходящийся с первым порядком. Для трехмерных задач на сетке, состоящей из шестигранников, возможно построение алго-
ритма, сходящегося только в слабом смысле, так как построить сопряженный объем, обеспечивающий выполнение условия аппроксимации потоков, явным образом не удается. Однако, если сетка состоит из тетраэдров, удается построить сильно сходящийся алгоритм. Алгоритм конструирования сопряженных объемов в этом случае аналогичен двумерному случаю.
Однако, в настоящей работе, алгоритм метода опорных операторов был использован только для задач фильтрации. Это связано с тем, что источники сил в задаче теории упругости имеют малые пространственные масштабы изменения, что приводит к плохообусловленным задачам. В связи с этим для задачи теории упругости был построен оригинальный алгоритм осреднения в ячейках для решения на сетках, более крупных, чем адаптированные к структуре среды. В основе алгоритма лежит требование равенства аппроксимации энергии. Как следствие, можно вывести слабую сходимость построенного алгоритма.
В настоящей работе в первой части рассмотрена постановка задачи Био. Кроме того, построена несвязанная математическая модель фильтрации в трещиноватых средах и эволюции их напряженно-деформированного состояния.
Во второй части рассмотрен метод опорных операторов как для двумерных скалярных эллиптических задач, так и для задач теории упругости при условии, что коэффициенты являются кусочно-постоянными. В качестве опорного оператора рассматривался оператор grad. Построена четкая теория этого метода для такого класса задач. Доказана сходимость в энергетической норме с тем же порядком, что и для постоянных коэффициентов.
Третья часть посвящена различным алгоритмам решения двумерных параболических уравнений на криволинейных сетках. Рассмотрена полностью неявная схема, а также регуляризованная схема. Выведено условие абсолютной устойчивости для регуляризованной схемы. Доказана абсолютная устойчивость и сходимость для полностью неявной схемы.
В четвертой части рассмотрены различные алгоритмы осреднения в ячейках для скалярных и векторных эллиптических задач: методы стандартного и узлового осреднения, а также метод Пабона, - как в двумерном, так и в трехмерном случае. Кроме изотропного случая, была рассмотрена полностью анизотропная среда для задачи теории упругости. Построение алгоритмов производилось как на Лебедевской, так и на произвольной другой сетке, состоящей из четырехугольников в двумерной задаче и шестигранников в трехмерном случае. Доказана слабая сходимость (в энерге-
тическом скалярном произведении) с первым порядком для построенных алгоритмов. В конце обсуждается вопрос сходимости построенных алгоритмов для полной несвязанной задачи пороу пру гости.
В последней пятой части представлены результаты тестовых расчетов. Во-первых, для параболических уравнений представлено сравнение регуля-ризованного алгоритма с полностью неявной схемой. Затем представлены результаты расчетов двумерных стационарных задач теории упругости по методам стандартного и узлового averaging'a на Лебедевских сетках. После этого показаны расчеты волн в трехмерных нестационарных задачах теории упругости по всем алгоритмам осреднения. Наконец, представлены результаты для трехмерной несвязанной задачи пороупругости по алгоритму, описанному в четвертой части.
Обобщение модели Био на случай трещиновато-пористой среды
Можно получить выражение для коэффициента рыхлости через дренажный модуль объемного сжатия пористой среды и модуль объемного сжатия скелета, Kd и Ks, соответственно. Так называемое дренажное деформирование происходит при постоянном давлении р = ро- Тогда пренебрегая постоянным членом — Сро % В определении (Tij и разлагая произвольную деформацию на всестороннее сжатие и сдвиг, получаем модуль сжатия для дренажных условий Kd = А + д.
Рассмотрим теперь деформирование насыщенного пористого тела, при котором тело подвергается всестороннему сжатию с постоянным давлением, равным давлению жидкости. С одной стороны, при таком деформировании оц = KskkSij = —p$ij, где Ка - модуль объемного сжатия скелета. С другой стороны, e = efcfc. Таким образом, Оц {Kd kk — СР)- Сравнивая полученные два выражения для компонент тензора напряжений, получаем, что коэффициент рыхлости Обобщение модели Био на случай трещиновато-пористой среды
Перейдем к рассмотрению трещиновато-пористой среды. Трещины в обычной среде можно описать как поры, но с существенной анизотропией. Присутствие же трещин в изначально пористой среде ведет к появлению двух пористостей и двух проницаемостей. Пусть mi - отношение объема пространства, занятого трещинами, к общему объему пространства, a vcii - отношение объема пор, к общему объему. В модели предполагается, что объем, занимаемый трещинами, существенно меньше объема пор. Таким образом, mi С тг Математическая модель фильтрации в трещиноватых средах, если не рассматривается влияние деформации на изменение трещиновато-пористого объема, имеет вид
Здесь Vi и Vi - скорости фильтрации в трещинах и порах, соответственно, которые определяются законом Дарси где pi и р2 давления жидкости в трещинах и порах соответственно.
Характерной особенностью неустановившегося движения жидкости в трещиноватых средах является наличие обмена флюидом между порами и трещинами. Именно этот процесс описывает последнее слагаемое в (1.11) и (1.12). Можно считать, что этот процесс квазистационарен, а также что скорость обмена определяется перепадом давлений жидкости в трещинах и порахр2— Pi, вязкостью жидкости и некоторыми геометрическими характеристиками среды. Из соображений анализа размерностей можно получить [Баренблатт и др., I960], что
Кроме того, в связи с тем, что гп\ -С тпз, уравнение (1.11) можно считать квазистационарным. Далее, как и для случая пористой среды, будем рассматривать слабосжимаемую жидкость а также предположим слабую зависимость обоих пористостей от давления т2 = mo2(l + /Зс2(р2 - Ро2))-Теперь рассмотрим процесс деформирования данной среды. Как и ранее, где
В силу того, что mi тг, мы приходим к такому же виду закона Гука, что и для обычной пористой среды. Окончательно, рассматриваемая модель Био для трещиновато-пористой среды имеет вид
Таким образом, рассматриваемая система состоит из параболической системы уравнений для давлений р\ и р2 и эллиптического векторного уравнения стационарной задачи теории упругости.
Очевидно, что система уравнений, описывающих эволюцию давлений в порах и трещинах, соответствует предположению о том, что давление в трещинах устанавливается гораздо быстрее, чем в порах. Это связано не только с тем, что гп\ С тг, но и с тем, что проницаемость к\ ; &2. Таким образом, система уравнений (1.16), (1.17) является жесткой. Упругое равновесие в системе устанавливается также мгновенно, поэтому рассматривается квазистатическая задача теории упругости. Еще одной особенностью задачи является существенное различие актуальных пространственных масштабов, поскольку размеры пористо-трещиноватой зоны много меньше размеров области, в которой формируется напряженное состояние. В силу этого является целесообразным применение методов осреднения в ячейках для решения задач теории упругости. Действительно, размеры ячейки при решении задач теории упругости заметно больше, чем размеры пористого коллектора. В силу уравнений Био градиент порового давления может рассматриваться как источник сил для задачи теории упругости. Так как размеры порового коллектора много меньше актуальных масштабов задачи, то возникает задача теории упругости с сосредоточенными источниками.
В следующей главе мы построим конечно-разностную схему на неортогональной сетке для эллиптического уравнения с разрывным коэффициентом. Эта схема является обобщением метода опорных операторов, который в [Самарский и др., 1996] был построен и подробно исследован лишь для постоянных коэффициентов. В той же главе мы также докажем сходимость в энергетической норме с первым порядком.
Метод опорных операторов в задаче теории упругости с разрывными коэффициентами
Далее нашей целью является нахождение величин даЬ и Ц. Их значения определяют способ вычисления интегралов в (2.1).
Из физических соображений ясно, что близость энергии и ее дискретного аналога, а также близость непрерывного потока через границу ячеек и его разностного выражения играют определяющую роль в проблеме точности данного метода. Близость потоков в соответствующей метрике влечет за собой сильную сходимость. С другой стороны, достаточным условием для слабой сходимости в смысле энергетического скалярного произведения является соответствующая аппроксимация интеграла энергии. Принципиально важно, что в двумерном случае при использовании метода опорных операторов на адаптированных сетках удается обеспечить выполнение обоих условий.
Так как разностное и непрерывное выражение для потока линейны, то для их равенства на некоторой линейной оболочке необходимо и достаточно потребовать их равенства для элементов базиса в этой линейной оболочке. В то же время, в силу того, что выражения для энергий квадратичны, необходимым и достаточным условием для их равенства является равенство всех попарных скалярных произведений. Линейная оболочка должна быть таковой, чтобы ее элементы с достаточной точностью аппроксимировали значения функции. Однако базис линейной оболочки не может содержать более, чем 4 вектора, исключая константу, (для двумерной задачи с сеткой из четырехугольников), так как иначе задача определения метрического тензора будет переопределена. Кроме того, метрический тензор должен содержать свободные параметры, которые могут быть использованы, чтобы обеспечить аппроксимацию потоков.
Пусть п - вектор нормали к линии разрыва коэффициента. Тогда п- это превалирующее направление изменения к в некоторой окрестности разрыва, то есть к = к(п г), где r-радиус-вектор. Пусть т - вектор, перпендикулярный тг, причем пара (т;п) положительно ориентирована. Далее предположим, что в области Г2 можно ввести криволинейные координаты, дважды непрерывно дифференцируемые по #, у, причем тип являются касательными к соответствующим координатным линиям. Таким образом, векторы тип можно определить в каждой точке области D, а линии разрыва локально могут быть аппроксимированы отрезками, что очень важно для дальнейших рассуждений.
Относительно точного решения будем предполагать, что нормальная к разрыву компонента потока &Ц и тангенциальная компонента градиента Jj достаточно гладкие, чтобы быть локально аппроксимированными константами в смысле некоторой нормы. То есть
Это означает, что локально в области Е р аппроксимируется функцией из ДЕ) = $pan{l,plyp2}. Функции из этой линейной оболочки представляют собой вариант естественного обобщения на случай разрывного коэффициента стандартных линейных функций из span{\,x,y}) которые обычно используются при постоянном к. В той области, где к непрерывно, ДЕ) переходит в линейное подпространство span{l7x,у}. Линейные функции могут быть использованы при построении разностных схем для выполнения как принципа тождества энергий, так и принципа тождества потоков [Самарский и др., 1996]. Они также используются при построении схем метода конечных элементов [Сьярле, 1980]. Однако в случае разрывного коэффициента данный базис не обеспечивает выполнения условия непрерывности потока на разрыве. Для этого базис должен быть перестроен. Именно это и сделано выше. Ниже изложена строгая математическая формулировка проблемы аппроксимации решений уравнений с разрывными коэффициентами.
Введем пространства Я1,1 (К, D), где К - некоторый тензор второго ран Пусть точное решение р Є HL,l(kI,D), где I - единичный тензор второго ранга. Рассмотрим для простоты случай, когда все линии разрыва представляют собой прямые. В этом случае в любой подобласти Т области D такой, что diam(T) = 0(h), существует аппроксимация точного решения функцией из L(T)t причем норма ошибки аппроксимации в Н1(Т) есть 0(h), а в Ьг(Г) - О(h2). А именно имеет место следующая теорема:
Теорема 2.1. (S.Moskow, V.Druskin et al., 1999; I.Babuska et al., 1994) Для любой функции p Є HL {kl,T) существует ее аппроксимация Ір Є ЦТ) и
Этот факт следует из стандартных свойств аппроксимации функций из Я1. В случае произвольной формы линий разрыва в силу их аппроксимируемости отрезками константа С будет еще зависеть от кривизны координатных линий тип.
Пусть Н - ячейка сетки М. Используя то, что в каждой ячейке Е точное решение можно хорошо аппроксимировать функцией из Ь(Е), потребуем для каждой ячейки Е сетки М точного совпадения непрерывной энергии /H k(Vp)2dV и разностной gabQaQb на функциях из L(E). Это даст нам возможность конкретизировать даЬ. Приводимые ниже три теоремы и их доказательства во многом схожи с результатами, опубликованными в [Самарский и др., 1996]. Пусть 1а - это вектор, по модулю равный длине ребра а и направленный вдоль положительного направления на этом ребре. Определим для функции р Є L(E) выражение интеграла энергии f- k{4p)2dV. Тогда справедлива следующая теорема. Здесь суммирование производится по всем узлам рассматриваемой ячейки, (1 а\ 1 ь)ф есть евклидово скалярное произведение векторов, сопряженных векторам 1а, 1ь с общей вершиной ф\ Бф есть некоторые площади, присоединенные к узлу ф и сумма которых равна площади Н ячейки Е; к есть значение коэффициента в ячейке Е. Доказательство: условие равенства энергий влечет равенство всех скалярных произведений: где рг,рі Є LCE), qla = Яа(рг)- Так как к в ячейке постоянен, то L(E) = span{l, х, у}. Тогда получаем, что даЬ должны удовлетворять соотношению где la,i - Ї-Я компонент вектора la. Частным решением этого уравнения является даЬ, определяемое (2.9). о Равенство энергий на некотором подпространстве, с помощью которого можно хорошо приблизить решение, влечет за собой слабую сходимость [S.Moskow, V.Druskin et al., 1999]. Однако, определив специальным образом площади Зф и Vi, можно аппроксимировать нормальные компоненты потока через границу области V . Из данной аппроксимации следует сильная сходимость в некоторой норме (теорема 2.4). Что касается площадей 5 , определим их следующим образом. Рассмотрим ячейку ABCD (рис. 2.1). Тогда, например, SA равняется половине площади треугольника ABD, SB равняется половине площади треугольника ABC и т.д. Данный выбор способа вычисления площадей Зф будет объяснен чуть позже. А сейчас получим некоторые свойства элементов даЬ. Выше мы говорили о выражении gabQaQb как об энергии и об gabQaPb как о скалярном произведении лишь формально. Теперь очевидно, что для определяемого формулой (2.9) даЬ выражение для энергии действительно положительно, и, тем самым, второе выражение есть скалярное произведение. Далее нам понадобятся еще несколько свойств даЬ, а также оператора G, определенного выше. Для этого сделаем более-менее разумные предположения относительно сетки: не будем рассматривать слишком "сплющенные" ячейки. А именно, пусть существуют константы са = 0(1) и С/ = 0(1), такие, что для любой пары ребер а и 6 с общей вершиной отношение их длин Ш удовлетворяет неравенству -f с\. А также, что все углы в ячейках были заключены между са и 7Г — са. Тогда справедлива следующая теорема
Метод опорных операторов в задаче теории фильтарции
В этой главе мы рассмотрим предложенный автором оригинальный метод аппроксимации для уравнений в неоднородных средах. Сначала будет рассмотрено скалярное эллиптическое уравнение в изотропной среде. Затем мы обобщим данный подход на систему уравнений теории упругости в изотропном и в полностью анизотропном случаях. В начале построение рассматриваемого алгоритма будет проведено для наглядности на Лебе-девских сетках, а затем мы сконструируем его аналог на обычных прямоугольных, а также произвольных криволинейных.
Один из возможных путей решения задач в неоднородных средах - это использование криволинейных сеток, адаптированных к структуре среды. Затем на них можно реализовывать конечно-разностные схемы, например, метод опорных операторов, или схемы метода конечных элементов. Однако при таком подходе возникают проблемы при обращении матриц в задачах с узкими зонами с контрастными параметрами.
Другой подход, который будет рассмотрен в этой главе, связан с использованием произвольных, удобных нам по какой-либо причине сеток, возможно неадаптированных к линиям разрывов. Однако в этом случае некоторые ячейки будут иметь достаточно сложную структуру. Для этих ячеек найдем эффективный тензор, который является полностью анизотропным в общем случае, даже если исходное уравнение было изотропным. Очень важно, чтобы эффективный тензор учитывал особенности исходной среды. Принцип тождества энергий является одним из фундаментальных для того, чтобы осредненная математическая модель имела свойства исходной среды. Он будет реализован для набора базисных векторов, учитывающих особенности среды. Число этих векторов определяют точность аппроксимации. В начале рассмотрим простейшее скалярное эллиптическое уравнение
Однако при переходе к произвольному расположению линий разрыва в многомерной задаче не существует очевидного обобщения этой формулы. Обычно, и даже в промышленных пакетах, задача upscaling a для многомерного случая решается алгоритмом, предложенным Кингом [King, 1989]. Он рассматривает прямоугольную ячейку с границами, параллельными осям координат, [Д ХІ+І] Х [yj;yj+i] и состоящую из четырех равных прямоугольников с постоянной проницаемостью (см. рис. 4.1). Для этой ячейки определяют диагональный эффективный тензор, находя эффективную проводимость некоторой электрической цепи, посредством использования аналогии между законом Дарси и законом Ома для электрической цепи.
Было показано [Максимов и др., 2002], что эффективный тензор по Кингу можно получить из решения конечно-разностными методами некоторых краевых задач на достаточно небольшом разностном шаблоне. А именно, для определения коэффициента кх необходимо рассмотреть граничные условия первого рода на вертикальных границах ячейки и нуль потока на горизонтальных границах. Неизвестные поместим в центры ячеек (см. рис. 4.1). Тогда, решив данную задачу, то есть определив pi, рг, рз и Р4 можно определить эффективный коэффициент кх приравняв два выражения для полного потока через ячейку:
Аналогичным образом необходимо поступить для нахождения ку. В случае произвольного расположения линий разрыва необходимо разбить ячейку на достаточно много подъячеек так, чтобы можно было считать, что в каждой подъячейке мы имеем структуру, используемую Кингом.
Однако при таком подходе можно определить лишь диагональные компоненты эффективного тензора. Попытки применить подобный подход для вычисления недиагональных компонент приводят к тому, что получаемый тензор оказывается несимметричным. Этот факт влечет за собой то, что соответствующая матрица системы линейных уравнений будет несимметрична, и стандартные методы обращения матрицы типа сопряженных градиентов неприменимы. Поэтому в дальнейшем нас будут интересовать лишь симметричные эффективные тензоры.
Вначале мы опишем так называемое стандартное осреднение [Moskow, Druskin et al., 1999]. Рассмотрим сначала простейший случай, когда линии разрыва параллельны оси у. Тогда Wx и меняются медленно в окрестности разрыва в то время, как Wy и 2 разрывны. Выражая быстро меняющиеся переменные через медленно меняющиеся и осредняя по ячейке простейшим интегрированием, мы получаем, что Таким образом,
Мы получили известный результат о том, что в двумерной задаче при рассматриваемом расположении линий разрывов эффективный тензор диагоналей, причем по одному направлению - это среднее гармоническое, а по другому - среднее арифметическое. В общем случае расположения линий разрывов применение данного метода, а также рассматриваемых в дальнейшем, возможно лишь на таких сетках, где в каждой ячейке имеют место только параллельные друг другу разрывы. В стандартном осреднении при таком ограничении эффективный тензор определяется как матрица квадратичной формы, которая в базисе, связанном с разрывом, представляется в виде Далее мы опишем алгоритм узлового осреднения в ячейках [Moskow, Druskin et ak, 1999] для указанной задачи. Нашей целью является построение конечно-разностной аппроксимации на Лебедевской сетке. Рассмотрим прямоугольную равномерную сетку: Лебедевская Р-сетка определяется как подмножество S с четной суммой i+j, а R-сетка состоит из оставшихся узлов S. Основным преимуществом Лебедевской сетки является то, что все компоненты потока определяются в одних точках по сравнению с обычной равномерной сеткой. Вследствие этого, оператор связи градиента с потоком, а также обратный к нему, являются локальными. Пока будем рассматривать, как в работе [Moskow, Druskin et ak, 1999], прямоугольные объемы, присоединенные к узлам и Р-сетки, и R-сетки. Со Рис 4.3. Расположение неизвестных в узлах Лебедевской сетки. Прямоугольные и ромбические ячейки. гласно терминологии, принятой в [В.П. Мясников, М.Ю. Заславский, А.Х. Пергамент, 2004], это - осреднение на прямоугольниках. Из рассуждений главы 2 следует, что точное решение в каждой ячейке Н R-сетки может быть аппроксимировано функцией из где тип- векторы, тангенциальный и нормальный к разрывам, соответственно. Мы определяем эффективный тензор Еу, используя принцип тождества энергий для функций из L(H). Это означает, что мы приравниваем непрерывную энергию, которая может быть вычислена аналитически, и разностную. А именно, так как непрерывный и разностный градиент р равен нулю, мы получаем я Здесь p,i - аппроксимация градиента р в R-ячейке Лебедевской сетки с помощью центральных разностей. Например для нумерации неизвестных, показанной на рис. 4.3, она имеет вид Таким образом, мы получаем три уравнения для трех компонент симметричного эффективного тензора модулей упругости. Теперь можно построить конечно-разностную схему для поставленной задачи
Трехмерная полностью анизотропная задача теории упругости
Математические модели в подземной гидродинамике, описывающие влияние потока флюида на напряженно-деформированное состояние окружающей среды, всегда представляли практический интерес. Например, при разработке гигантских углеводородных месторождений, характерных для России, возможно влияние процессов фильтрации на сейсмическую активность региона. При этом потоки флюида приводят к возникновению отсутствующих в геостатике касательных напряжений [Турунтаев и др., 1994].
Впервые подобная модель была предложена Терцаги [Terzaghi, 1936]. Он ввел эффективный тензор напряжений, зависящий от деформации скелета и смещений флюида. Однако применение соотношений Терцаги на практике представляет большие сложности. Реально могут быть наблюдаемы деформации порового объема как целого, но не деформации скелета и смещения флюида по отдельности. Сейчас процессы фильтрации в насыщенных пористых средах, сопровождаемые изменением напряженно-деформированного состояния, обычно рассматриваются в рамках квазистатической модели Био [Biot, 1962; Rice and Cleary, 1976]. В зависимости от того, рассматривается ли обратное влияние изменения напряженно-деформированного состояния на движение флюида, модель называется связанной или несвязанной.
В работах [Авербух и др., 1995; Гнедин и др., 1996] в рамках несвязанной модели Био были получены пространственные распределения давления фильтрующегося флюида, которые затем использовались для определения полей деформации и напряжений. В отличие от [Desrocher, 1994], где определялась только величина относительного изменения объема, в указанных выше работах были вычислены все компоненты тензора напряжений. В работе [А.В. Колдоба и др., 1999] была рассмотрена двумерная связанная задача Био, т.е. совместное решение уравнение фильтрации и уравнений теории упругости.
Однако в моделях такого типа пористая среда рассматривается как совокупность непроницаемых зерен, разделенных порами. В действительности же для всех естественных пластов характерна развитая в той или иной степени трещиноватость. Здесь возможны два подхода. В первом из них рассматривается более или менее регулярная система трещин. Однако, даже предположив, что мы сможем решить задачу с достаточно общим их распределением, все равно, находясь на поверхности Земли, невозможно достоверно судить о конфигурации этой системы. Здесь ситуация аналогична введению модели пористой среды. Даже если бы мы смогли проинтегрировать уравнения движения вязкой жидкости в порах, все равно их расположение в среде нам неизвестно. Так были введены осредненные величины, такие как пористость, давление, проницаемость, скорость фильтрации, а все основные законы формулировались в терминах этих величин. Именно по этому пути пошли сторонники второго подхода. Так была введена модель фильтрации в трещиноватых средах [Баренблатт и др., 1960].
В связи с вышесказанным представляет интерес построение единой модели фильтрации в трещиноватых средах и эволюции напряженно-деформированного состояния среды.
При этом увеличение давления приводит к появлению сдвиговых деформаций, которые могут стать причиной оползания пород на границе разлома. Такой процесс может спровоцировать сейсмическую активность среды вплоть до техногенного землетрясения.
Для описания фильтрации в трещиноватых средах используется модель с двойной пористостью и двойной проницаемостью [Баренблатт и др-, 1960]. Насыщенная среда представляет собой структуру, составленную из пористых блоков и трещин. Если макроскопическое описание пористой среды есть результат усреднения по масштабам большим размеров пор, то для трещиноватых сред вводятся характеристики, являющиеся результатом усреднения по масштабам большим сравнительно с размерами блоков. Отличие развиваемой здесь схемы от обычной схемы фильтрации в пористой среде состоит во введении в каждой точке пространства двух давлений жидкости - давления жидкости в порах и давления жидкости в трещинах -и учета обмена жидкостью между трещинами и порами. При определенных предположениях получается выражение для интенсивности этого обмена. Выводятся основные уравнения фильтрации жидкости в трещиноватой породе и более общие уравнения фильтрации жидкости в пористой среде с двойной пористостью.
Таким образом, уравнения фильтрации в трещиновато-пористой среде представляют собой уравнения с двойной пористостью и двойной проницаемостью. Причем, поскольку проницаемость в трещинах очень велика, давление в трещинах устанавливается гораздо быстрее, чем в порах, так что для его описания может быть использовано квазистатическое уравнение. Кроме того, можно предположить, что объем трещин пренебрежительно мал, и, как следствие, можно пренебречь величиной флюидосодержания в трещинах. В результате получается модель, отличающаяся от модели Био для пористых сред уравнением для давления в трещинах при условии, что обобщение закона Гука в модели Био сохраняет свое выражение. В дальнейшем мы используем несвязанный вариант данной модели, предложенный Николаевским [1996].
Хорошо известно, что для земной коры характерна структурная рас-слоенность. Эта расслоенность проявляется в чередовании зон высоких и низких скоростей, областей повышенных и пониженных напряжений. Существует множество примеров предельных случаев, когда фильтрация происходит в достаточно узких областях - разломах.
Необходимость учета геологической структуры среды, которая состоит из ряда слоев сложной геометрии, резко отличающихся по своим свойствам, привела к развитию численных методов решения разностных задач на криволинейных сетках, адаптированных к структуре среды. Построение конечномерной аппроксимации задачи на криволинейных сетках может быть осуществлено как проекционными методами, например, методом конечных элементов, так и методом конечных разностей. Если использовать метод конечных элементов с треугольными элементами, то в двумерном случае могут быть получены пятиточечные схемы для ортогональных сеток или семиточечные для неортогональных сеток при соответствующей ориентации прямоугольников. При соответствующем способе аппроксимации источников в такой задаче может быть доказана сходимость алгоритма в метрике Н1. В настоящей работе рассмотрен для данной задачи метод опорных операторов. Особенностью этого метода является использование при определении разностных операторов дивергенции и градиента с помощью разностного аналога известного интегрального тождества
