Н. Г. Чурбанова Москва 2011 г МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ МОСКОВСКИЙ ФИЗИКО-ТЕХНИЧЕСКИЙ ИНСТИТУТ (Государственный Университет) Факультет Управления и Прикладной Математики Д И П Л О М Н А Я Р А Б О Т А на тему: “ Математическое моделирование течений двухфазной слабосжимаемой жидкости в неоднородной пористой среде ” Студента гр. 571 Гичяна М.С. Научный руководитель: кандидат физико-математических наук, Н. Г. Чурбанова Москва - 2011 г. Содержание Введение……………………………………………………………………...…3 §1. Физическая постановка задачи …...…………………………………5 1. Общий вид модели …………….. ...............................................................5 2. Капиллярное давление и фазовые проницаемости в двухфазном случае………………. .......................................................................................7 §2. Математическая постановка ….....……………………………...…10 1 . Использование явных схем для моделирования процесса двухфазной фильтрации ……………………………………………………..................10 2.Математическая постановка двухфазной задачи...............................12 2.1 Система уравнений .....................................................................12 2.2 Разрыв среды ...............................................................................13 §3. Алгоритм решения…………………………………….………..14 1. Алгоритм решения в случае распространения загрязнения в однородной среде (двумерная область) …................................................14 2.Метод Ньютона...…................................................................................15 3. Алгоритм решения в случае неоднородной среды (двумерная область) .......................................................................................................17 2.1 Интерфейсные условия …..........................................................18 2.2 Разбиение расчетной области на группы узлов …..................19 §5. Результаты расчетов...................................................................25 1 . Просачивание пятна загрязнения в однородной среде.......................25 2 . Просачивание пятна загрязнения в среде со слабопроницаемой линзой ………………………………………………………………..….….28 Заключение ……………………………………………………………………33 Список литературы …………………………………………………......…..34 Введение. Теория фильтрации, изучающая законы движения жидкостей, газов и их смесей в пористой или трещиноватой среде, имеет обширное практическое применение. Фильтрационные расчеты играют исключительно важную роль при разработке технологий добычи нефти и газа, при проектировании, постройке и эксплуатации гидротехнических и мелиоративных сооружений, в горном деле, в решении экологических проблем. Не теряя своего значения на протяжении многих десятилетий, задачи фильтрации в настоящее время становятся еще более актуальными. В частности, теорию фильтрации можно применять для моделирования процессов, связанных с загрязнением почвы. Загрязнение почв нефтепродуктами или более плотными веществами относится к числу наиболее опасных, поскольку оно принципиально изменяет свойства почв, а очистка от них очень сильно затруднена. Загрязняющее вещество попадает в почву при различных обстоятельствах: при разведке и добыче нефти, при авариях на нефтепроводах, при авариях речных и морских нефтеналивных судов. Различные углеводороды попадают в почву на нефтебазах, бензозаправках и т.п. Последствия для почв, вызванные загрязнением, можно без преувеличения назвать чрезвычайными. Загрязняющее вещество обволакивает почвенные частицы, почва не смачивается водой, гибнет микрофлора, растения не получают должного питания. Наконец, частицы почвы слипаются, а само вещество постепенно переходит в иное состояние, ее фракции становятся более окисленными, затвердевают, и при высоких уровнях загрязнения почва напоминает асфальтоподобную массу. Бороться с таким явлением очень трудно. При малых уровнях загрязнения помогает внесение удобрений, стимулирующих развитие микрофлоры и растений. В результате нефть частично минерализуется, некоторые ее фрагменты входят в состав гуминовых веществ и почва восстанавливается. Но при больших дозах и длительных сроках загрязнения в почве происходят необратимые изменения. Тогда наиболее загрязненные слои приходится просто удалять. Настоящая работа посвящена проблеме, связанной с загрязнением почвы. Рассматривается двумерное течение жидкости сквозь пористую среду в вертикальном слое под воздействием силы тяжести. К загрязняющим веществам относятся минеральное топливо, растворители, очищающие средства (такие вещества называются NAPL – от английского Non-Aqueous Phase Liquids) [3]. В зависимости от того, как соотносится плотность вещества с плотностью воды, NAPL разделяются на легкие и плотные (первые называются Light NAPL, или LNAPL, их плотность меньше плотности воды; вторые – Dense NAPL, или DNAPL, соответственно их плотность больше плотности воды). Бензин, например, относится к первому типу, а тетрахлороэтилен – ко второму. Эти жидкости не смешиваются с водой и с газом, поэтому при моделировании их течения в почве говорят о многофазной фильтрации. В данной работе рассматриваются двухфазная фильтрация слабосжимаемых жидкостей . Первая фаза это вода, которой заполнен в начальный момент резервуар, вторая фаза – более плотная жидкость, которую ниже мы будем называть DNAPL. В работе было решено несколько тестовых задач. Двумерная задача о просачивани плотной жидкости через слабопроницаемую линзу с различными граничными условиями и мощностями источников DNAPL . Во всех задачах предполагается, что среда гетерогенная, то есть состоит из слоев с разными пористостью и проницаемостью. Так как на границе между разными материалами существует разрыв свойств среды, а именно абсолютной проницаемости и пористости, то на этой границе должны быть поставлены специальные интерфейсные условия для обеспечения физически корректного решения. Поскольку процесс просачивания не занимает длительного времени и происходит в относительно малом пространстве под воздействием только силы тяжести, при моделировании нужно рассматривать влияние капиллярных эффектов. Для решения задач подобного типа существует несколько подходов в зависимости от типа уравнений, образующих математическую модель. Традиционные подходы, такие как например IMPES метод, считают фильтрационные жидкости несжимаемыми и система уравнений для моделирования фильтрации в этом случае содержит эллиптическое уравнение для давления. В данной работе жидкости считаются слабосжимаемыми. В этом случае для описания процессов фильтрации обычно используется система уравнений параболического типа. Однако этот тип уравнения имеет достаточно жесткие ограничения на шаг по времени. В работе совершается переход к гиперболической системе уравнений, аналогично работам [5] и [6] , имеющей более мягкие ограничения на временной шаг. Для численной реализации используется явная трехслойная конечноразностная схема. Основной целью работы является разработка алгоритмов решения задач просачивания в почву со слоистой структурой (разной проницаемости и пористости) и исследование временной эволюции области загрязнения почвы. §1 . Физическая постановка задачи. 1. Модель в общем виде. Подземные пространства, в которых протекают гидродинамические процессы, представляют собой, как правило, пористые и трещиновато-пористые тела, характеризующиеся сложной системой пор, каналов, трещин, размеры которых малы по сравнению с характерными размерами среды, и по которым может происходить течение жидкости и (или) газа. Течение через такие пористые тела, при котором сила трения флюида (жидкости, газа) о скелет играет определяющую роль, называется фильтрацией. В работе проводится моделирование многофазной фильтрации несжимаемых несмешивающихся жидкостей. Под многофазной фильтрацией понимают совместное течение двух или более жидкостей через пористую среду. Количественной характеристикой пористости служит отношение объема пор к общему объему: (где m - коэффициент пористости, V п ­- объем пор, V – общий объем данного элемента среды). Фаза – это часть системы, которая является гомогенной (все ее компоненты имеют одинаковые физические и химические свойства). Физические или химические свойства фаз многофазной (гетерогенной) системы различны. Главными характеристиками движения многофазной системы является насыщенность и скорость фильтрации каждой из фаз. Насыщенностью Si порового пространства i -ой фазой называется доля объема пор, занятая этой фазой в элементарном объеме: . Очевидно, что . (1.1) Таким образом, в n -фазной системе имеется (n -1) независимая насыщенность. Поэтому при моделировании двухфазной фильтрации рассматривается лишь одна насыщенность. Центральную роль в теории фильтрации занимает понятие скорости фильтрации. Скорость фильтрации i -ой фазы в некоторой точке есть вектор такой, что объем данной фазы, протекающей за единицу времени через ориентированный элемент площади , расположенный в этой точке, включающий как поровое пространство, так и скелет, равен . Другими словами, проекция на некоторое направление равна отношению объемного расхода данной фазы к площадке, перпендикулярной к указанному направлению: . Предполагается, что справедлив обобщенный закон фильтрации Дарси [4], который для i -ой фазы при учете силы тяжести можно записать в виде: (1.2) где k – одна из феноменологичеческих характеристик пористой среды, называемая абсолютной проницаемостью, определяемая по данным о фильтрации однородной жидкости и не зависящая от свойств жидкости; μi – коэффициент динамической вязкости фаз; ki (Si ) – относительные фазовые проницаемости, определяемые экспериментально; Pi – давление в фазах; ρ i – плотность фаз; – вектор ускорения свободного падения; Si – насыщенность данной фазы. Для изменения со временем насыщенности какой-то фазы при фильтрации можно получить уравнение неразрывности, аналогичное уравнению неразрывности при обычном течении жидкости или газа (вне пористой среды). Рассмотрим некоторый объем Ω пористого тела. Масса флюида в этом объеме есть Масса флюида изменяется за счет потока флюида через границу ∂Ω, так что Считаем, что скелет неподвижен, тогда Получаем и в силу произвольности объема Ω приходим к уравнению неразрывности Если пористость m данной среды постоянна, то получаем такие уравнения для насыщенности фаз (1.3) Что приводит к классической модели фильтрации многофазной слабосжимаемой жидкости : где - плотность флюида, - давление, - скорость фильтрации, Ki и bi - соответственно коэффициенты проницаемости грунта и сжимаемости жидкости 2. Капиллярное давление и относительные фазовые проницаемости в двухфазном случае. Для задач, связанных с моделированием месторождений, например, при моделировании добычи нефти путём вытеснения водой, где время измеряется годами, а расстояния – километрами, влияние капиллярных сил на распределение давления незначительно и им можно пренебречь. Однако, при малых размерах области фильтрации и малых скоростях фильтрации капиллярные силы могут превзойти внешний перепад давления и их необходимо учитывать [5]. Капиллярные эффекты обусловлены межмолекулярными взаимодействиями двух различных фаз. Эти силы приводят к появлению угла смачивания на границе раздела двух фаз и к разрыву давления на этой границе. Разность фазовых давлений есть так называемое капиллярное давление : (1.4) где P 2 – давление несмачивающей фазы, а P 1 – давление смачивающей фазы. Жидкость с краевым углом θ < 90° называется смачивающей жидкостью относительно твердой фазы, а жидкость с краевым углом больше 90° – несмачивающей жидкостью (см. рис. 1). Рис. 1. Определение понятий смачивающей и несмачивающей фаз Относительные фазовые проницаемости ki , присутствующие в формуле (1.2), также обусловлены капиллярностью. Считается, что течение каждой из фаз происходит по своей системе связанных каналов, причем более смачивающая фаза занимает более узкие каналы. Несвязанные части каждой из фаз неподвижны. Из этого ясно, что относительные проницаемости зависят от насыщенностей Si . Величина капиллярного давления зависит от степени внедрения вытесняющей жидкости в пористую среду, поэтому капиллярное давление можно считать функцией от насыщенности смачивающей фазы: где – насыщенность смачивающей фазы (в данном случае - воды). Разность давлений двух фаз в одиночном капилляре радиуса r есть где σ – коэффициент поверхностного натяжения; θ – краевой угол смачивания. По аналогии, с учетом того, что средний размер пор есть , капиллярную разность давлений записывают в виде [5]: (1.5) Функция J (S ) называется функцией Леверетта. Она определяется структурой порового пространства и для каждого пористого тела, вообще говоря, своя. Но для группы пород сходного строения она одинакова. Достаточно хорошие приближения вида зависимости капиллярного давления от насыщенности дают модели Брукса – Кори (Brooks and Corey) и ван Генухтена (van Genuchten). Функция Брукса – Кори имеет вид , (1.6) где введена эффективная насыщенность , (1.7) где Swr – остаточная насыщенность смачивающей фазы в данной пористой среде, λ – показатель распределения размеров пор этой среды. В этой модели каждому пористому телу приписывается характеристическое значение P ­ d , называемое пороговым давлением. Pd – это минимальное капиллярное давление, необходимое для того, чтобы несмачивающая жидкая фаза начала вытеснять смачивающую фазу в соответствующей пористой среде. Зависимость капиллярного давления от насыщенности в модели ван Генухтена дается уравнением , где , то есть параметрами в этой модели являются α и n . При сравнении с уравнением Брукса – Кори видно, что α обратно пропорционально пороговому давлению P ­ d . Распределение размеров пор учитывается в параметре n . Большим значениям n соответствует более равномерное распределение, а малым (близким к 1) значениям – менее равномерное. Можно вывести уравнения, связывающие параметры этих двух моделей, но здесь они приводиться не будут. В данной работе используется модель Брукса – Кори. Относительные фазовые проницаемости также определяются в соответствии с моделью Брукса – Кори. Их зависимость от насыщенности разная для смачивающей и несмачивающей фаз и имеет вид (1.8) где Se определяется по формуле (1.7). §2 . Математическая постановка . 1.Использование явных схем для моделирования процесса двухфазной фильтрации. Рассмотрим классическую модель фильтрации однофазной жидкости: (2.1) При решении многих задач математической физики существуют минимальные размеры, меньше которых дальнейшая детализация описания уже не имеет смысла. В задачах фильтрации таким размером является величина l , порядок которой составляет несколько десятков диаметров зёрен породы. Учитывая этот минимальный размер l , в работах [7-8] по аналогии с кинетическими схемами [5] была предложена модель с модифицированным уравнением неразрывности: (2.2) где c – величина порядка скорости звука во флюиде. В отличие от исходной модели, модель с таким уравнением неразрывности позволяет использовать явные схемы с центрально-разностной аппроксимацией конвективного члена . В работах [7-8] показано, что для одномерной задачи, имеющей точное решение, полученное численное решение практически с ним совпадает . Заметим, что использование центральных разностей, имеющих второй порядок аппроксимации, дает преимущества при моделировании процесса фильтрации, особенно в геологических пластах сложной структуры. Вместе с тем, также как и для исходной модели, явные схемы для решения системы приводят к жёсткому ограничению ( )на допустимый шаг по времени. Представим разностную аппроксимацию производной по времени в виде: Опираясь на это представление, можно заменить уравнение неразрывности на гиперболическое уравнение (2.3) При шаге пространственного разбиения порядка h минимальный масштаб времени, меньше которого нет смысла проводить дальнейшую детализацию, имеет порядок (2.5) Здесь использовался тот факт, что скорость звука во флюиде c значительно превышает скорость фильтрации |u|. Таким образом, дополнительный член со второй производной в уравнении (2.3) по порядку величины значительно меньше, чем : (2.6) где L – характерный размер задачи. Тогда разностное уравнение для (2.3) примет вид: Обратим внимание, что в случае, когда справедливо равенство схема перейдёт в схему Дюффорта-Франкела [9], обладающую абсолютной устойчивостью. К сожалению, точность расчётов по схеме Дюффорта-Франкела в многих случаях является неудовлетворительной. Однако в рассматриваемом случае параметр специально выбирается из соображений, чтобы дополнительный член со второй производной по времени в уравнении вносил минимальные изменения в решение уравнения. С учётом (2.5) из (2.6) следует условие устойчивости для трёхслойной схемы решения уравнения: Таким образом, с учётом минимальных масштабов по пространству и времени исходная модель превратилась в модель: (2.7) которая позволяет использовать явные схемы с центрально-разностной аппроксимацией конвективных членов с более мягким условием устойчивости. 2. Математическая постановка двухфазной задачи. 2.1 Система уравнений Опираясь на результаты параграфов 1.2 и 2.1, выпишем систему уравнений для модели двухфазной фильтрации: (2.1) (2.2) (2.3) (2.4) (2.5) (2.6) (2.7) (2.8) (индекс α обозначает фазу: w – вода, n - NAPL от английского Non-Aqueous Phase Liquid) Где m – это пористость, Sα - насыщенность фазы, kα (Sw ) – относительная фазовая проницаемость, g – ускорение свободного падения, qα - источник жидкости, pc (Sw ) – капиллярное давление,Se - это эффективная насыщенность, Swr - остаточная насыщенность водной фазы, Pd - пороговое давление, λ – индикатор пористости для данной среды Относительные фазовые проницаемости определяются в соответствии с моделью Брукса – Кори. Их зависимости от насыщенности имеют вид (2.8) Рис. 2. Капиллярное давление по Бруксу-Кори с параметрами λ и Pd Рис. 3. Относительные проницаемости по Бруксу-Кори 2.2 Разрыв среды При рассмотрении процессов в неоднородной среде (то есть, среде, состоящей из блоков с различными свойствами – абсолютной проницаемостью, пористостью и т. д.) возникает проблема разрыва некоторых параметров на границе раздела фаз. В этом случае при расчетах необходимо выполнение двух условий на интерфейсе: непрерывность потока через границу раздела и непрерывность капиллярного давления. Из последнего условия вытекает разрыв насыщенности на интерфейсе (см. рис.4): Рис. 4. Непрерывность капиллярного давления и разрыв насыщенности на интерфейсе. . §3 . Алгоритм решения . 1. Алгоритм решения в случае распространения загрязнения в однородной среде(двумерная область). Выпишем разностное уравнение, получаемое для трехслойной явной схемы, для уравнения (2.1), в двумерном случае: (3.1) А также разностное уравнение для нахождения компонент скорости фильтрации: (3.2) (3.3) Общий вид расчетной области представлен на рис.1 Рис.5 Тогда алгоритм расчетов будет состоять из следующих этапов по каждой фазе: 1. Вычисление по формулам (3.2) и (3.3) компонент скорости фильтрации. 2. Вычисление из уравнения (3.1) на новом слое по времени. 3. Вычисление насыщенности DNAPL и давления водной фазы на новом слое по времени в каждом узле расчётной сетки из следующей системы уравнений: (3.4) Система (3.4) получена умножением уравнений состояний (2.3) с обеих сторон на и применением формулы для зависимости давлений фаз и капиллярного давления (2.5) и может быть решена методом Ньютона. 4. Вычисление , и из (2.4), (2.5) и уравнения состояния (2.3) соответственно. 2.Метод Ньютона Пусть требуется решить систему уравнений (4.1) где — заданные, вообще говоря, нелинейные (среди них могут быть и линейные) вещественные функции n вещественных переменных Обозначив , , данную систему (2.1) можно записать одним уравнением (4.2) Пусть ( ) — некоторая последовательность невырожденных вещественных n x n-матриц. Тогда, очевидно, последовательность задач , k = 0,1,2,... имеет те же решения, что и исходное уравнение (2.1а), и для приближенного нахождения этих решений можно формально записать итерационный процесс , k = 0,1,2,... (4.3) имеющий вид метода простых итераций (4.2.1b) при . В случае - это действительно МПИ с линейной сходимостью последовательности ( ) Если же различны при разных k , то формула (4.3) определяет большое семейство итерационных методов с матричными параметрами . Рассмотрим некоторые из методов этого семейства. Положим , где — матрица Якоби вектор-функции F(x). Подставив это в (4.3), получаем явную формулу метода Ньютона , (4.4) Эту формулу, требующую обращения матриц на каждой итерации, можно переписать в неявном виде: . (4.5) Применение (4.5) предполагает при каждом k = 0,1,2,... решение линейной алгебраической системы относительно векторной поправки , а затем прибавление этой поправки к текущему приближению для получения следующего: . К решению таких линейных систем можно привлекать самые разные методы как прямые, так и итерационные в зависимости от размерности n решаемой задачи и специфики матриц Якоби (например, можно учитывать их симметричность, разреженность и т.п.). Сравнивая (4.5) с формальным разложением F(x) в ряд Тейлора , видим, что последовательность ( ) в методе Ньютона получается в результате подмены при каждом k=0,1,2,... нелинейного уравнения F(x) = 0 или, что то же (при достаточной гладкости F(x) ), уравнения линейным уравнением т. е. с пошаговой линеаризацией. 3. Алгоритм решения в случае неоднородной среды(двумерная область). Общий вид расчетной области представлен на рис.2: Рис.6 В верхней части области находиться источник. А также в области присутствует среда с меньшей проницаемостью. (3.5) Как уже говорилось выше на границе раздела сред должны быть реализованы интерфейсные условия. В данной работе они реализованы на верхней границе раздела сред, а на боковой и нижней границах параметры сред усреднены следующим образом: (3.6) Расчетная сетка строиться таким образом, чтобы граница раздела сред проходила через узлы счета это дает возможность корректно построить решение на «интерфейсе», а также дает возможность не усреднять проницаемость K на других границах раздела сред. 3.1 Интерфейсные условия Рис. 7. Расчетная ячейка на границе раздела сред Если расчетная область состоит из двух подобластей с разными свойствами, то, как было сказано выше, капиллярное давление считается непрерывным на границе раздела, а водонасыщенность является разрывной величиной. Для моделирования этой ситуации предлагается следующий алгоритм. Обозначим соседние подобласти G 1 и G 2 , причем в G 2 пороговое давление больше, чем в G 1 . Тогда водонасыщенность в подобласти G 2 на интерфейсе рассчитывается по следующей формуле: (3.7) где – функция, обратная функции Брукса – Кори в области G 2 , выражающая зависимость насыщенности от капиллярного давления; и – значения водонасыщенности на интерфейсе соответственно со стороны областей G 1 и G 2 ; – водонасыщенность, при которой достигается пороговое давление области G 2 , то есть она определяется из условия (см. рис. 8). Таким образом, чтобы DNAPL попал в область G 2 , его насыщенность в области G 1 на интерфейсе должна быть больше, чем . (3.8) = (3.9) Рис. 8. Непрерывность капиллярного давления и разрыв насыщенности на интерфейсе. 3.2 Разбиение расчетной области на группы узлов Так как разностная схема является скрытой пятиточечной по каждой оси , то для корректного вычисления скоростей фильтрации можно разбить расчетную область на группы узлов как показано на рис. 9 Рис.9 Разбиение расчетной области на группы узлов (расчетные узлы находятся в серединах ячеек) 1-Область расчетные узлы, которой и смежные с ними узлы, полностью лежат в более проницаемой среде - . 2-Область расчетные узлы, которой лежат в более проницаемой среде , но имеют смежные узлы принадлежащие границе раздела двух сред - . 3-Область расчетные узлы, которой принадлежат верхней границе раздела сред . 4-Область расчетные узлы, которой принадлежат левой границе раздела сред- . 5-Область расчетные узлы, которой принадлежат нижней границе раздела сред- . 6-Область расчетные узлы, которой принадлежат правой границе раздела сред- . 7-Расчетный узел, находящийся на левом верхнем угле раздела сред- . 8- Расчетный узел, находящийся на правом верхнем угле раздела сред- . 9- Расчетный узел, находящийся на левом нижнем угле раздела сред- . 10- Расчетный узел, находящийся на правом нижнем угле раздела сред- . 11- Область расчетные узлы, которой лежат в менее проницаемой среде, но имеют смежные узлы принадлежащие границе раздела двух сред- . 12- Область расчетные узлы, которой и смежные с ними узлы, полностью лежат в менее проницаемой среде- . 1) Счет в первой и двенадцатой области Так как все точки разностной схемы для расчета для точек первой и двенадцатой области принадлежат одной среде, то в них расчет ведется аналогично счету в однородной области (см. §3 пункт 1). 2) Счет во второй и одиннадцатых областях Так как расчетные узлы этих областей имеют смежные лежащие на границе раздела сред, то в них нельзя напрямую как в §3 пункт 1, корректно рассчитать скорости фильтрации хотя бы по одному из направлений. Поэтому предлагается в точках этих областей перейти к уравнению неразрывности в интегральной форме. Приведем алгоритм расчета для второй области(для одиннадцатой расчет аналогичен, заменяются только параметры среды ) Уравнение неразрывности примет вид: (3.10) Тогда разностное уравнение для расчетного узла примет вид: (3.11) Скорости фильтрации в полуцелых точках вычисляются по формулам: (3.12) Плотности в полуцелых относительная проницаемость в полуцелых точках усредняются следующим образом: (3.13) Далее, выражая , подставляем получившиеся значения в уравнения (3.4) и приходим, как и в §3 пункте 1, к системе уравнений, которую можно решать методом Ньютона. Далее алгоритм аналогичен описанному в §3 пункте 1. 3) Счет в третьей, седьмой и восьмой областях В расчетных узлах данных областей ставятся интерфейсные условия, следовательно в данных узлах присутствует разрыв давлений хотя бы одной из фаз и для корректного расчета, так же предлагается перейти к записи уравнения неразрывности в интегральной форме. Приведем алгоритм расчета для третьей области (для седьмой и восьмой алгоритм аналогичен) Рассмотрим расчетный узел интерфейса (Рис. 10). Уравнение неразрывности в этом случае примет вид: (3.14) Где и - значение функции над и под интерфейсом Рис.10 Расчетный узел интерфейса Так как = const до просачивания через интерфейс. Следовательно, пока разностное уравнение будет иметь вид: (3.15) В данном уравнении - это значение функции над интерфейсом, скорости фильтрации находятся аналогично скоростям из (3.12) с учетом того какой среде принадлежат точки по которым они рассчитываются. Для плотностей и относительных проницаемостей применяются формулы аналогичные (3.13).Член уравнения потоку вектора через границы ячейки узла счета. Когда станет меньше чем произойдет переход DNAPL через интерфейс и разностное уравнение измениться следующим образом: (3.16) Где и -значении функции над и под интерфейсом соответственно Из формул (3.7) и (3.9) получаем: (3.17) Раскладывая данную функцию в ряд Тейлора до слагаемого с первой производной и подставляя в (3.16), а так же используя тот факт, что обе жидкости слабосжимаемы, приходим к разностному уравнению: (3.18) Где: (3.19) Далее, подставляя в уравнения (3.4), переходим к решению системы уравнений методом Ньютона. Затем аналогично §3 пункту 1находим значения всех функций над интерфейсом. Так как для обеспечения интерфейсного условия достаточно разрывности давления одной из фаз, полагаем непрерывным на интерфейсе. Затем из уравнений (3.7) и (3.17) находим насыщенности под интерфейсом, а остальные неизвестные величины под интерфейсом находятся из уравнений (2.5), (2.6) и (2.3). 4) Счет в четвертой, пятой, шестой, девятой и десятой областях Расчетные узлы данных областей лежат на границе раздела сред и свойства сред в данных узлах усредняются, но при этом присутствует разрыв в проницаемости сред k, который тоже можно усреднить, но в данной работе предлагается перейти к интегральной записи уравнения непрерывности, так же как и в предыдущих пунктах. Алгоритм построение разностного уравнения для нахождения такой же, как и в предыдущем пункте с той только разницей что все величины кроме проницаемости k, считаются непрерывными и на границе раздела усредняются по формулам: Так как рассмотрены все точки расчетной области и следовательно построен алгоритм расчета значений насыщенности фаз по насыщенностям на двух предыдущих слоях по времени. Алгоритм устойчиво считает при , что согласуется с теорией, описанной в §2 пункте 1. §5. Результаты расчетов 1. Просачивание пятна загрязнения в однородной среде Рассматривается двумерная задача о просачивании DNAPL в резервуар, заполненный водой. Параметры среды: Среда Единицы измерения m 0.4 – k 6,64·10–1 1 м2 Swr 0.09 – λ 3.86 – Pd 755 Па μ1 = 0.001 кг/(м·с); μ2 = 0.0009 кг/(м·с); ρ1 = 1000 кг/м3 ; ρ2 = 1460 кг/м3 ; В середине верхней границе поставлен источник q = 5·10–8 кг/с. а) начальные условия . б) граничные условия (боковые границы): (верхняя граница): (нижняя граница): Размер расчетной области - 30x50 расчетных узлов(по y и x соответственно) Размер узла 0.1м x 0.1м На рисунках 11 а, б и в показано распределение насыщенности DNAPL в различные моменты времени. Рис.11a Распределение насыщенности DNAPL в момент времени 720сек. Рис.11 б Распределение насыщенности DNAPL в момент времени 2160 сек. Рис.11 в Распределение насыщенности DNAPL в момент времени 3600 сек. S y На рис. 12 показан график зависимости (т.е в центральном срезе) в момент времени 3600 сек. Рис.12 График зависимости 2. Просачивание пятна загрязнения в среде со слабопроницаемой линзой Рассматривается двумерная задача о просачивании DNAPL в резервуар, заполненный водой. В центре резервуара находиться слабопроницаемая линза. Среда Линза Единицы измерения m 0.4 0.4 – k 6,64·10–1 1 7,15·10–12 м2 Swr 0.09 0.12 – λ 3.86 2.49 – Pd 755 2060 Па μ1 = 0.001 кг/(м·с); μ2 = 0.0009 кг/(м·с); ρ1 = 1000 кг/м3 ; ρ2 = 1460 кг/м3 ; В середине верхней границе поставлен источник q = 5·10–7 кг/с. а) начальные условия . б1) граничные условия (условия не протекания) (боковые границы): (верхняя граница): (нижняя граница): Размер расчетной области - 50x71 расчетных узлов(по y и x соответственно) Размер узла 0.1м x 0.1м На рисунках 13 а, б и в показано распределение насыщенности DNAPL в различные моменты времени. Рис.13a Распределение насыщенности DNAPL в момент времени 960сек. Рис.13б Распределение насыщенности DNAPL в момент времени 2880сек. Рис.13в Распределение насыщенности DNAPL в момент времени 4800сек. На рис. 14 показан график зависимости (т.е в центральном срезе) в момент времени 4800 сек. S y Рис.14 График зависимости б2) граничные условия( q = 3·10–5 кг/с) (боковые границы): (верхняя граница): (нижняя граница): Размер расчетной области - 20x35 расчетных узлов(по y и x соответственно) Размер узла 0.1м x 0.1м На рисунке 15 показано распределение насыщенности DNAPL в момент времени 1900 сек. (Линза и источник не симметричны относительно центра области) Рис.15 Распределение насыщенности DNAPL в момент времени 1900сек. На рис. 16 показан график зависимости (т.е в центральном срезе) в момент времени 1900 сек. Рис.16 График зависимости Заключение. В дипломе построена модель вертикального просачивания двухфазных слабосжимаемых жидкостей под воздействием силы тяжести сквозь неоднородную среду (состоящую из областей с различными проницаемостью, пористостью и т. д.). Создан алгоритм расчета, использующий специальные условия на границе раздела двух сред. Были проведены расчеты задачи о просачивании плотной жидкости сквозь слабопроницаемую линзу в резервуаре, наполненном водой. В дальнейшем возможно усовершенствования модели путем проведения расчетов с использованием неравномерных сеток на многопроцессорных вычислительных комплексах. Список литературы. 1. Chetverushkin B.N. High-performance computing: Fundamental problems in industrial application / In Parallel, distributed and grid computing for engineering. Stirlingshire: Saxe-Coburg Publications, 2009, p.369–388. 2. Басниев К.С., Кочина И.Н., Максимов В.М. Подземная гидромеханика. - М.: Недра, 1993. 3. Четверушкин Б.Н. Кинетически-согласованные схемы в газовой динамике. - М.: Изд-во МГУ, 1999. 4. Четверушкин Б.Н. К вопросу об ограничении снизу на масштабы в механике сплошной среды // Время, хаос, математические проблемы, вып. 4, Изд-во МГУ, 2009, с.75-96. 5. Трапезникова М.А., Белоцерковская М.С., Четверушкин Б.Н. Аналог кинетически-согласованных схем для моделирования задачи фильтрации // Математическое моделирование, 2002, т.14, №10, с.69-76. 6. Четверушкин Б.Н., Морозов Д.Н., Трапезникова М.А., Чурбанова Н.Г., Шильников Е.В. Об одной явной схеме для решения задач фильтрации // Математическое моделирование, 2010, т.22, №4, с.99-109. 7. Самарский А.А. Теория разностных схем. - М.: Наука, 1977. 8. Азиз Х., Сеттари Э. Математическое моделирование пластовых систем. - М.: Недра, 1982.