автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Вычислительные алгоритмы для задач однофазной и двухфазной фильтрации на основе схемы КАБАРЕ

Московский государственный университет им. М.В. Ломоносова Факультет вычислительной математики н кибернетики

На правах рукописи

Канаев Антон Андреевич

Вычислительные алгоритмы для задач однофазной и двухфазной фильтрации на основе схемы КАБАРЕ

005008467

(05.13.18- Математическое моделирование, численные методы и комплексы программ)

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

1 9 ЯНВ 2012

Москва—2011

005008467

Работа выполнена в Институте проблем безопасного развития атомной энергетики Российской академии наук (ИБРАЭ РАН)

Научный руководитель доктор физ.-мат. наук, профессор

Головизнин Василий Михайлович

Официальные оппоненты доктор физ.-мат. наук, профессор

Фаворский Антон Павлович

доктор физ.-мат. наук, профессор Цыпкин Георгий Геннадиевич

Ведущая организация: Институт вычислительной математики РАН

Защита состоится « » СрСЯрСИЛ, 201Д^ода в/Р.ЗОчасов в аудитории №685. на заседании диссертационного совета _Д 501.001.43 в Московском государственном университете имени М.В. Ломоносова по адресу: 119991, ГСП-1, Москва, Ленинские горы, МГУ имени М.В. Ломоносова, 2-й учебный корпус, факультет ВМК.

С диссертацией можно ознакомиться в библиотеке факультета ВМК МГУ имени М.В. Ломоносова.

Автореферат разослан «_»_2011г.

Ученый секретарь

диссертационного совета ¿Э

д.ф.-м.н, профессор / Захаров Е.В.

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ1

Актуальность темы

Тема диссертации связана с проблемой безопасности атомной энергетики, в частности, с безопасным захоронением долгоживущих радиоактивных отходов. Рассматриваемые в настоящее время схемы замыкания ядерного топливного цикла не предполагают стопроцентного выжигания продуктов ядерного деления, и из всех возможностей их утилизации наиболее предпочтительной является технология долговременного (на сотни тысяч лет) подземного захоронения. Безопасность таких захоронений зависит от интенсивности процессов переноса растворенных солей в геологических формациях мест их расположения. Большую актуальность приобретает задача разработки математических моделей миграции радиоактивных растворов на временах, сравнимых с геологическими. Разработке математических моделей и новых вычислительных алгоритмов для расчета влагоперсноса в сильно неоднородных трещиновато-пористых средах и посвящена настоящая диссертация.

История вопроса

Математические модели, описывающие динамику многокомпонентных флюидов в пористых и трещиноватых геологических породах с учетом реальных уравнений состояния и термодинамических процессов в настоящее время хорошо известны и широко применяются при расчетах геотермальных источников и задач нефте- и газодобычи (Баренблатг Г.И., Коновалов А.Н, Пергамент А.Х., Цыпкин Г.Г., К. Pruess, J.P. Gwo, Селяков В.И.). Несколько другой класс моделей используется при решении задач фильтрации промышленных и экологических загрязнений через зону неполного влагонасыщения в грунтовые воды.

Общим для этих классов задач является то, что результаты, полученные с точностью до десяти процентов, с практической точки зрения являются вполне приемлемыми. Совершенно иные требования предъявляются к результатам анализа безопасности захоронения радиоактивных отходов. Здесь требуются иные постановки задач и новые методы их решения.

В задачах подземной гидродинамики часть уравнений, описывающих перенос концентрации, являются гиперболическими. Это определяет возможность существования особенностей решения в виде сильных и слабых разрывов и супердиффузионной асимптотики в скорости распространения «хвостов» загрязнения. Расчеты задач просачивания в сильно неоднородных трещиновато-пористых средах представляют собой серьезную вычислительную проблему, поскольку наличие у большинства известных алгоритмов аппроксимационной вязкости (в рассматриваемом случае «аппроксимационных капиллярных сил») существенно искажает характер решения в экстремальных случаях. Естественным решением, в этой ситуации, представляется выбор численной схемы, относящейся к т.н. алгоритмам «высокой разрешающей способности»(А. Harten, S. Osher, В. van Leer, P.D. Lax, X.-d. Liu, C.-W. Shu, T. Chan, Фаворский А.П., Бакирова М.И., Тишкин В.Ф., Вязников K.B., Карпов В.Я.). Достаточно полный обзор работ этого направления содержится в монографии (Куликовский А.Г., Погорелов Н.В., Семенов А.Ю., 2001) Альтернативным подходом можно считать использование схемы КАБАРЕ.

Схема КАБАРЕ для простейшего одномерного линейного уравнения переноса была предложена и подробно исследована в работах Головизнина В.М. и Самарского A.A. в 1998 году. Позже выяснилось, что в западпой литературе она известна как схема Айзерлиса (Upwind Leapfrog) поскольку является представителем семейства разностных

'Работа выполнена при частичном финансировании по грантам РФФИ №0608-01501-а и № 06-0 l-00819-a.

3

схем, исследованных эти автором на устойчивость в 1986 году. В дальнейшем, эта схема претерпела ряд эволюционных скачков. Значительный вклад в ее развитие внесли Головизнин В.М., Карабасов С.А., Кобринский И.М. К наиболее важным изменениям исходной трехслойной схемы КАБАРЕ (Upwind Leapfrog) можно отнести ее представление в двухслойном виде, что было осуществлено введением дополнительных неизвестных, т.н. «консервативных» переменных. Затем она была дополнена алгоритмом коррекции потоков, базирующемся на прямом использовании принципа максимума. Это сделало возможным ее обобщение на более содержательные и сложные по сравнению с простейшим уравнением переноса, случаи.

Для задач подземной гидродинамики схема КАБАРЕ бьша впервые использована в работе (Головизнин В.М., Карабасов СЛ.,1998) для решения системы уравнений Лаверетга-Бакли в двумерной задаче о скважине. Работа носила чисто методический характер - в ней было проведено сравнение четырех возможных способов обобщения схемы КАБАРЕ, три из которых основывались на трехслойной версии алгоритма. Результат, полученный по одной из этих схем, оказался на то время удовлетворительным. Дальнейшее развитие схемы КАБАРЕ привело к необходимости критического пересмотра полученных ранее результатов с целью их качественного улучшения на новой методической основе.

Целью настоящей работы является.дальнейшее развитие схемы КАБАРЕ для решения одномерных нелинейных законов сохранения с выпуклыми и невыпуклыми функциями потоков; разработка на базе модифицированной схемы нового эффективного вычислительного алгоритма для численного решения одномерных задач однофазной и двухфазной фильтрации в зонах неполного влагонасьпцения; разработка новых математических моделей двумерной однофазной и двухфазной фильтрации в геологических формациях с сильно неоднородными фильтрационными свойствами.

Научная новизна

На основе сформулированного в работе принципа минимума парциальной локальной вариации консервативных переменных получен новый вид численного потока для одномерного квазилинейного гиперболического уравнения, который можно рассматривать как нетривиальную модификацию численного потока Годунова для схемы КАБАРЕ. Показано, что новые численные потоки позволяют эффективно решать не только задачи с выпуклыми функциями потоков, но и задачи с невыпуклыми потоками, такими как потоки в уравнении Лаверегга-Бакли. Решения получаются монотонными, ударные волны размазываются на одну-две расчетные ячейки. Показано, что схема КАБАРЕ с новыми потоками, при включении в нее дополнительного регулируемого диссипативного механизма, дает возможность получать истинные энтропийные решения для случаев, в которых другие алгоритмы высокой разрешающей способности, такие как TVD схемы не способны вычленять решения с максимальной энтропией, что приводит к расцеплению комплекса из ударной волны и волны разрежения и их контакту через зону постоянного течения.

Для одномерной задачи о просачивании влаги в зоне аэрации разработаны два новых эффективных вычислительных алгоритма на основе схемы КАБАРЕ, один из которых основывается на однофазном представлении и учете сил реакции трещиновато-пористой среды (модель просачивания в вакууме), а второй - на численном решении по схеме КАБАРЕ одномерной двухфазной системы уравнений Лаверетга-Бакли с учетом силы тяжести и поверхностного натяжения. В первой модели не учитывается эффект вытеснения воздуха, во второй влияние этого фактора учтено в полной мере. Показано, что при соотношении подвижностей воздуха и воды е =0.01 влияние вытесняемого воздуха на процесс просачивания незначительно и им можно пренебречь. При этих условиях модель однофазного просачивания с учетом сил реакции породы в вычислительном плане оказывается более эффективной.

Схема КАБАРЕ обобщена на двумерную квазилинейную систему двухфазных уравнений Лаверетга-Бакли, что позволило существенно уменьшить паразитные вычислительные эффекты, связанные с наличием псевдокапиллярных сил, имеющих аппроксимационную природу.

Разработана новая математическая модель двумерной однофазной фильтрации в зоне неполного влагонасыщения с учетом сил реакции геологической среды. Предложено «гидростатическое» приближение этой модели для численного исследования динамических процессов просачивания влаги по перколяционным решеткам. Использование гидростатического приближения позволяет отказаться от решения сильно вырожденных разностных уравнений эллиптического типа па давление, что приводит к существенной экономии вычислительных ресурсов.

Достоверность результатов работы обеспечивается четкостью постановок задач и строгостью математических доказательств, многочисленными тестовыми и модельными расчетами на сгущающихся расчетных сетках, сравнением полученных результатов с аналитическими решениями и результатами других авторов.

Теоретическая и практическая значимость

В предыдущих версиях схемы КАБАРЕ для одномерного квазилинейного скалярного гиперболического уравнения для определения направления экстраполяции потоковых переменных использовалось направление скорости звука, что приводило к неэнтропийным решениям в случаях сложных невыпуклых функций потока с множественными экстремумами. Отказ от использоваши малодостоверной информации о направлении прихода характеристики и формулировка задачи о вычислении новых потоковых переменных с использованием принципа минимума парциальной локальной вариации консервативных переменных привела к повышению надежности расчетов с невыпуклыми функциями потоков.

Новые вычислительные алгоритмы для решения задач однофазной и двухфазной фильтрации в зоне неполного влагонасыщения дают возможность получения несглаженных решений в случае геологических сред с сильно контрастными коэффициентами пористости и проницаемости, что позволяет адекватно моделировать пористые среды с фрактальными сетями трещин.

На защиту выносятся

Новый вид численных потоков дня схемы КАБАРЕ. Исследования пределов применимости новой модификации методики КАБАРЕ для задач с невыпуклой функцией потоков.

Новые вычислительные алгоритмы для одномерной однофазной и двухфазной фильтрации в зоне неполного влагонасыщения на основе схемы КАБАРЕ.

Новые математические модели двумерной однофазной и двухфазной фильтрации в зоне аэрации.

Публикации и апробация работы.

Основные результаты были представлены на следующих конференциях и семинарах:

- XVII Всероссийская конференция "Теоретические основы и конструирование численных алгоритмов решения задач математической физики", посвященная памяти К.И. Бабенко, Абрау-Дюрсо, 2006

- X Школа Молодых Учёных ИБРАЭ РАН, Москва, 22-23 апреля 2009

- XI Школа Молодых Учёных ИБРАЭ РАН «Безопасность и риски в энергетике», Москва, 22-23 апреля 2010

- Международная научно-техническая конференция «Суперкомпьютерные технологии: разработка, программирование, применение» СКТ-2010,

Дивноморское, Гелевджикский район, Краснодарский край, лечебно-оздоровительный комплекс «Голубая даль», 27 сентября-2 октября 2010г. - XII Школа Молодых Учёных ИБРАЭ РАН «Безопасность и риски в энергетике» Москва, 28-29 апреля 2011

Основные результаты диссертации опубликованы в 3 статьях, в том числе в материалах 1 международной конференции. 1 статья размещена в издании, рекомендованном ВАК.

Структура и объем работы

Диссертация состоит из введения, трех глав основного текста, заключения и списка литературы.

КРАТКОЕ СОДЕРЖАНИЕ РАБОТЫ

Во введении рассмотрен основной круг проблем, связанных с темой диссертации, приведен краткий обзор литературы, сформулирована цель работы, защищаемые положения и научная новизна.

В первой главе диссертации описывается новый алгоритм численного решения одномерного скалярного квазилинейного гиперболического уравнения

Bu^dF(u) dt "

+ —^ = 0 (1)

с достаточно произвольной функцией конвективных потоков .

К уравнению такого вида в простейшем случае сводится задача о просачивании влаги в зоне неполного влагонасьпцения.

Решению такого типа задач посвящены работы отечественных ученых Кружкова С.Н., Ладыженской O.A., Фаворского А.П., Головизнина В.М., Тишкина В.Ф., Олейник O.A., Васильева В.И., Василевского Ю.В. и др.

Теория построения «энтропийных» решений одномерных квазилинейных гиперболических уравнений была развита Кружковым С.Н. во второй половине двадцатого века. Известен для таких уравнений и метод нахождения автомодельных решений задачи Римана. Однако общие методы численного решения систем гиперболических уравнений, базирующиеся на точных или приближенных методах решения задачи Римана о распаде произвольного разрыва надежно работают только для выпуклых функций Сходимость схем с низкой диссипацией (схемы высокой разрешающей способности) к энтропийному решению критическим образом зависит от деталей алгоритма подсеточного восполнения и для произвольной невыпуклой функции, является скорее исключением, чем правилом.

Новый метод, предложенный в диссертации, базируется на схеме КАБАРЕ, оперирующей двумя типами переменных. Первые относятся к центрам расчетных ячеек и называются «консервативными», вторые - «потоковые», относятся к расчетным узлам. Таким образом, на текущем временном слое считаются заданным как «консервативные», так и «потоковые» переменные.

Расчетный цикл в схеме КАБАРЕ можно разложить на фазы «предиктор-корректор». На первой фазе по явной схеме вычисляются промежуточные консервативные переменные на полуслое по времени. Затем линейной экстраполяцией значений потоковых переменных через промежуточную консервативную находятся потоковые

переменные на новом временном слое. Завершает цикл фаза «корректор», на которой происходит вычисление новых значений консервативных переменных.

Схема КАБАРЕ устойчива при числах Куранта от нуля до единицы, обладает вторым порядком аппроксимации на неравномерных пространственно-временных расчетных сетках, бездиссипативна во всей зоне устойчивости и является точной при двух числах Куранта - единице и одной второй. Кроме того, схема КАБАРЕ определена на максимально компактном вычислительном шаблоне - все арифметические операции в ней осуществляются только с величинами в пределах одной пространственно-временной расчетной ячейки.

Недостатком линейной схемы КАБАРЕ является ее немонотоннсть. Для придания численному решению свойства монотонности новые значения потоковых переменных повергаются нелинейной коррекции, базирующейся на принципе максимума. Вычисленные по формуле линейной экстраполяции их величины не должны выходить за пределы интервала между максимальным и минимальным значениями искомой функции на нижней грани пространственно-временной расчетной этейки.

Ранее, в публикациях Головизнина В.М. и Кобринского И.М. была предпринята пытка обобщения схемы КАБАРЕ на случай нелинейных функций потоков. В простейшем случае она сводилась к тому, что направление экстраполяции потоковых переменных ( из «левого нижнего угла ячейки в правый верхний» или «из правого нижнего - в левый верхний») выбиралось в зависимости от направления скорости звука в центре расчетной ячейки. Для выпуклых функций такой алгоритм приводил к хорошим результатам. Однако, для невыпуклого потока Лаверетга-Бакли дела обстояли не столь благополучно -при распаде разрыва на волну разрежения и ударную волну параметры последней не выходили на энтропийное решение. Для улучшения ситуации было предложено использовать в каждой расчетной ячейке анализатор эволюционности разрыва

В работах Головизнина В.М. и Кобринского И.М. описана реализация этой идеи, которая привела к удовлетворительным результатам при решении уравнения Лаверетта-Бакли, однако в более сложных случаях проблемы все еще оставались. Это и послужило основанием для разработки нового универсального подхода, представленного в первой главе.

Поскольку скорость звука для невыпуклых функций перестает быть надежным ориентиром для выбора направления диагональной экстраполяции, то последнюю предлагается проводить в обоих направлениях, и результаты корректировать в соответствии с принципом максимума. Для каждого расчетного узла на новом временном слое, таким образом, будет получено два разных значения потоковой переменной - одно со стороны левой прилегающей ячейки, которое будем: обозначать как и"1, другое со стороны правой ячейки - и"^. Основной вопрос заключается в выборе на основе имеющейся информации адекватного численного потока Нас

При построении разностных схем высокой разрешающей способности, таких как ТУЛ и ТУВ решение о выборе численных потоков выбираются из анализа их влияния на полную вариацию решения. В первом случае численные потоки выбираются из условия уменьшения полной вариации, во втором - из условия ограниченности ее роста. Этот хорошо зарекомендовавший себя подход положен в основу выбора численного потока и в настоящей работе.

Рассматривается локальная вариация консервативных переменных на промежуточном временном слое, прилегающих к данному узлу расчетной сетки

Предлагается выбирать такие значения и".*1 потоков!,IX переменных из допустимого интервала [«"'>""г]* чтобы соответствующий поток 17(и"*')приводил к быстрейшему уменьшению локальной вариации, а если это невозможно, то к ее минимальному росту.

Сформулированный подход к определению численного потока назван «принципом парциального минимума локальных вариаций».

Далее показывается, что для выпуклых (впуклых) функций К(и) искомое значение численного потока Я0С(а,А) определяется выражением:

Нетрудно видеть, что полученное выражение представляет собой модификацию известного потока Годунова.

На примере многочисленных тестовых расчетов показано, что модифицированные потоки Годунова приводят к энтропийным решениям как при выпуклых, так и при невыпуклых функциях Г (и). Рассмотрены случаи потоков Лаверетга-Бакли,

^"(м) = 8т(и), ^(и)=|(м-1)2-1|идр.

Вторая глава диссертации посвящена некоторым математическим вопросам, возникающим при численном решении уравнения просачивания влаги в сильно неоднородных геологических формациях. Эту задачу можно формулировать как двухфазную, когда первой фазой является вода, а второй - вытесняемый ей воздух. Такая постановка приводит к системе из двух уравнений Лаверетта-Бакли, описывающих баланс каждой фазы. В случае, когда второй фазой является нефть, система уравнений Лаверетга-Бакли достаточно хорошо изучена и широко применяется на практике. Плотности и подвижности воды и нефти относительно близки, поэтому трудностей принципиального характера при численном решении этих уравнений не возникает, либо они достаточно легко преодолеваются. Иное дело, когда второй фазой является воздух или газ, отличающийся по этим параметрам от воды на несколько порядков.

Рассмотрим задачу об одномерном просачивании в грунт при резких разрывах коэффициента проницаемости (сильно контрастных геологических средах). Для этого рассмотрена двухфазная одномерная модель Лаверетта-Бакли без учета сил поверхностного натяжения:

где а = к"1; Ь = и\

<Г'

/и-г^ + Лу (№„) = <); т——~ + Лу(Д11) = 0

&

(2)

и соответствующие законы Дарси:

V/.

2

Р.-*];

Здесь s(x,t) - функция влагонасыщенноста, т (х) - пористость среды, к(х) -абсолютная проницаемость, р^.р,, - плотности воды и воздуха соответственно, g -гравитационное ускорение, p{x,t) - давление, /,(*),/а(з)в - относительные проницаемости воды и воздуха, dv,,da - подвижности воды и воздуха. Для воды и воздуха

djdw=t= ю-:.

Показано что для одномерного случая задача сводится к решению одного квазилинейного уравнения вида (1) без каких либо дополнительных связей:

ds 5F(x,s) п т— + —^—i = 0;

8t дх

F(x,s)=*^.L(s).[l-G(S)] + Q(t)-G(sy, (4)

G(i) [/.W+-/-W]' e4'

По схеме, предложенной в первой главе, проведены: численные расчеты модельной задачи при разных соотношениях подвижностей воздуха и воды. В ячейках, где влагонасьпценносгь приближается к единице, возникает проблема резкого уменьшения шага по времени для получения устойчивых результатов, и, соответственно, увеличения расчетного времени.

Естественно рассмотреть предельный случай, когда плотность второй субстанции стремится к нулю, а ее подвижность - к бесконечности, что соответствует просачиванию влаги сквозь пористое тело в вакууме. Если ограшгеиться одномерным случаем и отвлечься от процессов капиллярного смачивания, то в такой идеализированной постановке проблема сводится к решению начально-краевой задачи для уравнения

8s 8F(x,s) „

—+—^ = 0 (5)

8t дх w

с необязательно выпуклой функцией F(x,s) при наложенной на решение ограничивающей связи

о s si (6)

Четвертый раздел второй главы посвящен исследованию применимости алгоритма КАБАРЕ, предложенного в первой главе, для функций потока, зависящей от искомой функции и пространственной координаты х без учета ограничивающей связи (6). На примерах решения тестовых задач показано, что схема КАБАРЕ приводит к энтропийным решениях! при наличии сильных разрывов в функции F(x,s).

В следующем разделе описывается модификация схемы КАБАРЕ, ориентированная на учет ограничения (6). Предполагается, что на текущем временном слое это ограничение выполнено. Делается шаг по времени в соответствии с алгоритмом КАБАРЕ и находятся консервативные величины концентраций на новом шаге по времени. Попутно вычисляются и новые значения потоков. Если все консервативные концентрации удовлетворяют ограничению (6), то решение на новом слое считается найденным. Если хотя бы в одной ячейке это ограничение нарушается, включается итерационная процедура коррекции вычисленных на первом этапе потоков. В результате возникает недоопределенная система уравнений на неизвестные вариации потоков, из которой

нельзя получить однозначное решение. Для нахождения единственного решения привлекается дополнительное условие - принцип минимума вариаций потоков. Выписывается квадратичный функционал от искомых вариаций, и решается задача его условной минимизации с учетом связей, возвращающих «зашкалившие» консервативные переменные в нужный диапазон. Это приводит к необходимости численного решения системы линейных алгебраических уравнений относительно неизвестных множителей Лагранжа со специфически вырожденной трехдиагональной матрицей, которая решается прямым методом Гаусса. Показано, что возникающие при решении вышеописанной задачи множители Лагранжа имеют физический смысл сил реакции механической системы на отклонение от заданного ограничения и могут бьпь интерпретированы как силы давления.

По вычисленным в результате первой итерации силам давления находятся вариации потоков, которые приводят к перераспределению консервативных значений концентраций. При этом зафиксированные ранее отклонения от связей будут устранены, но могут возникнуть отклонения в других ячейках, в которых ранее их не было. Если это так, то проводится следующая итерация. И так до тех пор, пока условие (6) не будет выполнено для всех расчетных ячеек. Работоспособность описанного алгоритма демонстрируется на ряде тестовых примеров.

Сравнение полученных результатов с полученными ранее с учетом влияния воздуха показало качественно схожие результаты.

Капиллярные силы часто играют в процессе просачивания влаги в пористую среду определяющую роль, поэтому их необходимо учитывать в соответствующих математических моделях. Уравнение баланса влаги с учетом сил поверхностного натяжения можно представить в виде модифицированного уравнения «конвекции -диффузии»

& д 81 дх ~ дх

< \ &

(7)

где пористость, для упрощения, полагается равной единице и у Способ учета

диффузионной правой части этого уравнения в вычислительный алгоритм схемы КАБАРЕ описан в последнем разделе второй главы диссертации.

Однопараметрическое семейство схем второго порядка аппроксимации для предиктора и заключительной части корректора записывается как:

И

т/2 А

т/2 И

где оператор £ определяется выражением:

¿•У

. 1

V (у г" \ ^МД'^111

(8)

(9)

(Ю)

а параметр со е [0,1].

Уравнение (8) естественным образом расщепляется на два;

т/2 А 1 '

5п+1/2 оя+1/2

М/2 £ ^н+1/2

т/2

Первое из этих уравнений является явным и легко разрешается. Для решения второго, при <з*0, следует использовать специальную форму алгоритма прогонки. Блок экстраполяции потоковых переменных и их нелинейной коррекции в алгоритме КАБАРЕ, описанный в предыдущих разделах диссертации, при учете диффузии остается неизменным.

Линейный анализ устойчивости описанного алгориша показывает, что допустимый шаг по времени должен определяться из условия:

__„ . , 1 к

Т ^ СИ-ПШИТ ,Т 1, Т„=---¡=-7—Г-;-=г,

' 1 J (13)

= _1 _#_

' 2 тах[а,-у(«;)/т,>1/2]

для явной схемы (ш = 0), и

тйС^- хг (14)

для неявной (о = I).

Известно, что функция у при приближении s к единице неограниченно возрастает. Условие (13) в этом случае оказывается излишне обременительным и в расчетах приходится применять неявную схему с параметром и = 1. При большом разбросе значений коэффициентов оператора (10) для численного решения неявной разностной схемы (12) необходимо использовать потоковую прогонку.

(И) (12)

В третьей главе диссертации результаты, полученные в первых двух главах, применяются для построения базового фрагмента двумерной модели просачивания влаги в сильно неоднородную пористую среду с трещинами, имеющими фрактальную природу. Описываются новые методы расчета просачивания влаги по перколяционным сетям и пористым матрицам с сильно неоднородными коэффициентами пористости и проницаемости.

Первый раздел третьей главы посвящен обобщению алгоритма КАБАРЕ для уравнений Лаверетга-Бакли, описанного во второй главе, на случай двух пространственных переменных и демонстрации специфических свойств этого алгоритма.

Рассматривается система уравнений двухфазной подземной гидродинамики без учета сжимаемости флюидов и породы, сил поверхностного натяжения и силы тяжести:

где л - концентрация воздуха (газа), ф (?) - пористость. Эта система приводится к известному виду

ф ^' 5+)=^ (Г' *) ■ &**(р)]=-(Л+/,),

об)

РгИ,

представляющему собой одно гиперболическое уравнение относительно концентрации газообразной фазы и одно эллиптическое уравнение относительно давления. Рассматривается разностная аппроксимация этих уравнений на расчетной сетке с квадратными ячейками. Как и в предыдущих главах, вводятся консервативные и потоковые переменные концентрации газовой фазы, аппроксимируется закон сохранения концентрации, алгоритм записывается в форме «предиктор-корректор», для вычисления потоковых переменных используется линейная экстраполяция, дополненная нелинейной коррекцией на основе принципа максимума.

Аппроксимация второго уравнения приводит к стандартному пятиточечному разностному уравнению относительно давления.

Новая система явно-неявных разностных уравнений аппроксимирует систему уравнений (16) со вторым порядком как по времени, так и по пространственным переменным. Выписывается условие устойчивости, характерное для явной схемы КАБАРЕ.

Все ранее известные обобщения схемы КАБАРЕ на уравнения газовой динамики и гидродинамики несжимаемой жидкости демонстрировали резкое уменьшения диссипативных явлений, обусловленных ошибками аппроксимации, по сравнению с другими, ранее известными схемами. Для демонстрации аналогичного свойства у схемы КАБАРЕ для задач подземной гидродинамики рассмотрена модельная задача о вытеснении жидкой фазы газовой в квадратной области при наличии непроницаемой трещины «5-6» (рис 1).

Рис.1

Рис. 2

В заданной области система уравнений (16) решается со следующими начальными и граничными условиями: на верхней границе области задано постоянное давление, равное единице, на нижней границе - постоянное давление, равное нулю, на двух боковых границах и на горизонтальной трещине - условие непротекания. В начальный момент концентрация газа во всей области равна нулю.

Относительная подвижность газа и жидкости полагается, как и в предыдущей главе, равной б = рг • ^/р, =0.01. Относительные проницаемости задаются как

На рис. 2 в изометрии изображено поле давления. Видно, в окрестности непроницаемой трещины поле давления испытывает разрыв. На поверхности поля давления светлым цветом обозначена концентрация газа. Фронт вытеснения, отодвигается от верхней границы, и начинается развитие его неустойчивости, приводящее к возникновению пальцеобразных струй (фингеринг). При наличии аппроксимационных эффектов, характерных для классических вычислительных алгоритмов, прогалины между струями быстро затягиваются, и длина границы вытеснения резко сокращается. Схема КАБАРЕ позволяет свести этот паразитный эффект к минимуму.

Второй раздел третьей главы посвящен обобщению однофазного приближения теории просачивания на многомерный случай.

Задача рассматривается в двумерной постановке. Предполагается, что ускорение свободного падения направлено по оси х. Уравнение, описывающее динамику влаги, в этом случае имеет вид:

(17)

ду1

ду

<Гу(х'У)-гП

& Эу]

где

к (*>>■.')=<*Лх'У)

8 Р*~

, ку(х,у,1) = -ау(х,у)

др ду'

Цх.у)

а -

Решение уравнения (17) должно удовлетворять ограшиениям:

0 £ .г (х, у, г) 21

(18)

При просачивании влаги по сети относительно крупных трещин влиянием поверхностного натяжения можно пренебречь и уравнение (17) упрощается:

Следует особо подчеркнуть, что основные вычислительные проблемы вызывают т.н. «сильно контрастные среды», когда пористость и проводимость представляют собой случайные разрывные функции с вариациями, стремящимися к бесконечности.

Обобщение схемы КАБАРЕ на двумерный однофазный случай также разбивается на несколько этапов.

На первом этапе вычисляются консервативные переменные влагонасьпценности на промежуточном временном слое по явной схеме первого порядка аппроксимации по

13

времени и второго по пространственным переменным. На втором этапе происходит вычисление предварительных значений потоковых неременных на гранях расчетных ячеек по формулам линейной экстраполяции для вертикальных граней и для горизонтальных.

Далее происходит нелинейная коррекция полученных экстраполяцией потоковых переменных, основанная на использовании принципа максимума с модификацией на ненулевую правую часть для каждого направления экстраполяции и, наконец, решение для каждой грани расчетной ячейки задачи о распаде разрыва с модифицированными потоками Годунова.

Консервативные переменные на новом временном слое вычисляются на третьем этапе алгоритма. Если все они удовлетворяют условию (18), то решения считается найденным и происходит переход на следующий временной слой. В противном случае вычисляются отклонения консервативных концентраций от допустимых границ (невязки) и решается вырожденное уравнение Пуассона относительно неизвестных давлений с вычисленными невязками в правых частях. Этот процесс повторяется до тех пор пока не будет выполнено условие (18).

При решении уравнения (19) с ограничениями (18) на перколяционных решетках проницаемость граней которых либо нуль, либо заданная фиксированная величина,' описанный выше алгоритм существенно упрощается за счет более простого алгоритма учета ограничений. При возникновении на пути потока влаги непреодолимого препятствия, выстраивается гидростатический столб, давление которого и создает необходимую для выполнения ограничений (18) реакцию. Соответствующее приближение в работе названо «гидростатическим». Использование гидростатического приближения на порядок уменьшает требования к вычислительным ресурсам и существенно упрощает вычислительный алгоритм.

Завершают главу примеры решения тестовых и модельных задач.

Рассмотривается задача о протечке влаги на перколяционной решетке размером (50,50) ячеек. Относительная проницаемость берется в виде квадратичной функции Хопфа, а пористость считается постоянной и равной единице. Исходное влагосодержание в области равно нулю. На левой границе единичного квадрата поддерживается постоянное влагосодержание = 0.4, на верх ей и нижней границах - нулевые потоки, на правой

границе - нулевое значение влагосодержания. Перколяционный порог, при превышении которого на двумерной решетке образуется бесконечный перколяционный кластер, отвечает перколяционному параметру т|' » 0.59. Рассматриваются различные реализации случайных полей при т! = 0.5, 0.55 и 0.6. Последний случай смоделирован также на более подробной решетке (200x200) при числе Куранта СРЬ = 0.45.

Динамика протечки влаги с левой границы области до правой изображена на рис. За-£

Рис. За. Влагосодержание, ИТ = 30000 Рис. ЗЬ. Давление, ЫТ = 30000

14

ОСНОВНЫЕ РЕЗУЛЬТАТЫ ДИССЕРТАЦИИ

- Разработан новый вид численных потоков для схемы КАБАРЕ. Проведены исследования пределов применимости новой модификации методики КАБАРЕ для задач с невыпуклой функцией потоков. Показано, что включение в алгоритм диссипэтора Паниковсхого с коэффициентом \ = 0.6 позволяет получать энтропийные решения для невыпуклых функций потока при высоких порядках касания производной в точках смены выпуклости. Установлено, что к такому же эффекту приводит сужение области возможных значений потоковых переменных, допускаемой принципом максимума.

- На основе схемы КАБАРЕ разработаны новые вычислительные алгоритмы для одномерной однофазной и двухфазной фильтрации в зоне неполного влагонасыщения. Задача однофазной фильтрации сформулирована как задача решения одномерного квазилинейного гиперболического уравнения при наличии ограничивающей связи.

- Разработан алгоритм учета такой связи в алгоритме КАБАРЕ с использованием вариационного принципа для возмущений потоков. Предложен способ учета капиллярных эффектов. Для двухфазной фильтр ции с учетом явления вытеснения воздуха получено аналитическое выражение для градиента давления через величину влагонасыщения, что привело к классической одномерной гиперболической задаче с немонотонной функцией потоков. Проведено сравнение численных результатов решения модельной задачи о просачивании по двум этим моделям, Показано, что при реальных значениях отношения подвижностей воздуха и воды физические результаты получаются близкими, в то время как вычислительные затраты в однофазной модели существенно меньше.

- Разработаны новые математические модели двумерной однофазной и двухфазной фильтрации в зоне аэрации. Для двумерной двухфазной модели Лаверетта-Бакли с учетом вытесняемого воздуха разработан вычислительный алгоритм на основе схемы КАБАРЕ. Вычислительные особенности этого алгоритма продемонстрированы на модельной задаче о напорном течении в пористой среде с непроницаемой трещиной.

- Представлена двумерная постановка задачи об однофазной фильтрации с учетом ограничений на решение. Разработан эффективный вычислительный алгоритм, являющийся обобщением одномерного алгоритма. Для убыстрения счета задач о динамике протечек на перколяционных решетках разработано т.н. «гидростатическое» приближение. Работоспособность и вычислительная эффективность нового подхода продемонстрирована на ряде модельных и тестовых задач.

ПУБЛИКАЦИИ АВТОРА ПО ТЕМЕ ДИССЕРТАЦИИ

- «Новый вычислительный алгоритм для математического моделирования просачивания влаги сквозь ненасыщенную трещиноватую геологическую среду с такой проницаемостью», Головизиин В.М., Семенов В.Н., Канаев A.A. и др. - М 2006. - 53 с. - (Препринт / ИБРАЭ РАН; NIBRAE-2006-07)

- «Новый вычислительный алгоритм для математического моделирования просачивания влаги сквозь ненасыщенную трещиноватую геологическую среду с низкой проницаемостью», Канаев A.A., Материалы Международной научно-

технической конференции «Суперкомпьютерные технологии: разработка, программирование, применение», 27 сентября - 2 октября 2010 г., с. «Принцип минимума парциальных локальных вариаций для определения конвективных потоков при численном решении одномерных нелинейных скалярных гиперболических уравнений», Головизнин В.М., Канаев A.A., Ж. вычисл. матем. и матем. физ., 51:5 (2011), 881-897

Введение.

Глава 1. Модификация конвективных потоков в схеме КАБАРЕ для одномерных нелинейных скалярных гиперболических уравнений.

1.1 Схема КАБАРЕ для простейшего линейного уравнения переноса.

1.2 Проблема переключения потоков в схеме КАБАРЕ.

1.3 Обобщение схемы КАБАРЕ на случай нелинейных потоков.

1.4 Частная задача Римана для уравнения с выпуклыми потоками.

1.5 Форма представления оператора Римана, не опирающаяся на свойство дифференцируемости функции потока.

1.6 Процедура согласования начальных значений консервативных и потоковых переменных и оператор переключения потоковых переменных

1.7 Невыпуклые функции потоков. Принцип минимума парциальной локальной вариации.

1.8 Одномерные квазилинейные уравнения с произвольными потоками. Примеры расчетов.

Глава 2. Одномерные модели двухфазной и однофазной фильтрации в зоне неполного влагонасыщения.

2.1 Постановка задачи о просачивании влаги в зоне аэрации.

2.2 Одномерное двухфазное просачивание. на основе схемы КАБАРЕ, и был сформулирован «принцип минимума парциальных локальных вариаций» для определения конвективных потоков при численном решении одномерных нелинейных скалярных гиперболических уравнений.

2.3 Просачивание влаги в пустую пористую среду.

2.4 Скалярный закон сохранения с функцией потоков, зависящей от координат.

2.5 Алгоритм учета удерживающих связей.

2.6 Задача о просачивании жидкости в пустую (заполненную вакуумом) пористую среду с учетом капиллярных сил.

Глава 3. Двумерные модели двухфазной и однофазной фильтрации в зоне неполного влагонасыщения.

3.1 Прецизионный явно-неявный алгоритм для решения системы двумерных уравнений Лаверетта-Бакли с использованием схемы КАБАРЕ

3.2 Тестовая задача о дифракции фронта вытеснения воды (нефти) газом на непроницаемой трещине.

3.3 Обобщение однофазного приближения теории просачивания на многомерный случай.

3.4 Протекание по двумерным перколяционным решеткам. Гидростатическое приближение.

Образование радиоактивных отходов — неотъемлемое свойство ядерных технологий, а безопасное обращение с ними — важная часть проблемы безопасности ядерной энергетики. Основная задача здесь состоит в том, чтобы исключить распространение в окружающей среде радиоактивных веществ, образующихся при работе атомных станций, предотвратить их воздействие на человека и природу при хранении.

В ближайшие десятилетия наиболее реальным технически достижимым методом, способным обеспечить изоляцию отработанного ядерного топлива и остеклованных высокоактивных отходов от окружающей среды в течение 10 тыс. или более лет, является удаление этих ядерных отходов в глубокозалегающие породы — гранит, глина, соль, туф.

Многие из нуклидов, которые содержатся в продуктах переработки ядерного топлива, имеют очень большой период полураспада, исчисляемый сотнями и тысячами лет, а для распада трансурановых нуклидов (актинидов) нужно время порядка миллиона лет. Понятно, что сегодня у людей нет опыта контроля над какими-либо техническими устройствами в течение столь длинного периода времени. Решить данную проблему можно с помощью надежного прогнозирования поведения радионуклидов в подземном хранилище. Один из подходов к решению данной задачи заключается в построении математических моделей, адекватно описывающих перенос радионуклидов в сильно неоднородных, неупорядоченных геологических структурах.

Математические модели, описывающие динамику многокомпонентных сред в пористых и трещиноватых геологических породах с учетом реальных уравнений состояния и термодинамических процессов в настоящее время хорошо известны и широко применяются при расчетах геотермальных источников и задач нефте- и газодобычи [1-6]. Несколько другой класс моделей используется при решении задач фильтрации промышленных и экологических загрязнений через зону неполного влагонасыщения в грунтовые воды [1-14]. Вычислительные методики, реализующие эти модели, также хорошо известны [1, 2, 15-26].

Общим для этих классов задач является то, что результаты, полученные с точностью до десяти процентов, с практической точки зрения являются вполне приемлемыми. Совершенно иные требования предъявляются к результатам анализа безопасности захоронения радиоактивных отходов. Здесь требуются иные постановки задач и новые методы их решения.

При анализе безопасности долговременного хранилища радионуклидов на полигоне штата Невада (Yucca Mountain) возникает задача моделирования протечек в водоненасыщенном пласте вулканического туфа, обладающего малой проницаемостью.

Зона неполного насыщения (зона аэрации). Основной вопрос, касающийся оценки проводящих свойств ненасыщенных трещиноватых пород — это вопрос о механизме распространения влаги в таких средах. Согласно классической капиллярной модели [27, 28], вода за счет капиллярного эффекта впитывается в твердую матрицу и распространяется за счет фильтрации по ней. При этом трещины являются препятствиями для движения воды на большие расстояния. Для сред, наиболее интересных с точки зрения задачи захоронения отходов (например, трещиноватый туф Yucca Mountain), по-видимому, более реалистична другая модель, согласно которой основным каналом распространения воды являются именно трещины, а капиллярное впитывание представляет сравнительно слабый эффект [2, 29]. В этом случае фильтрационный поток оказывается крайне неоднородным и нестационарным, наблюдаются наличие предпочтительных путей распространения. Режим распространения воды и переноса примесей определяется статистическими свойствами сети трещин.

К настоящему времени проведено немного экспериментальных исследований фильтрации воды и транспорта примесей в трещиноватых ненасыщенных породах. В [30, 31] описана серия полевых экспериментов с измерениями структуры фильтрационного потока и транспорта примесей через трещиноватую породу, а также внутренней структуры трещин. Наблюдения [30, 31] показали пространственную и временную нестабильность потока, сильное каналирование, когда большая часть потока (70-100%) проходит через малую часть доступного сечения трещин (15-20%). При этом «активные» пути движения воды постоянно меняются в зависимости от циклов смачивания/осушения, химического взаимодействия потока со стенками, отложений на стенках растворенных в воде материалов. Все это значительно усложняет картину фильтрации в ненасыщенной зоне по сравнению с насыщенной.

В [32] описаны результаты пневматических испытаний, проведенных в туннеле на Yucca Mountain, Nevada, USA, который является предположительным местом захоронения отходов. В блоке с размерами 15x20x15 м было просверлено в разных направлениях около 30 скважин длиной 5-10 м. В каждой скважине имелся уплотненный участок, через который подавался воздух с постоянным расходом. Одновременно с подачей воздуха измерялось давление в самой нагнетаемой скважине и во всех остальных. Процедура нагнетания и измерения давления повторялась последовательно со всеми скважинами. Отклик давления быстро, в течение нескольких минут устанавливался на стационарном значении. Эти данные позволили оценить проницаемость вблизи каждой скважины (усредненную по длине уплотненной области), а также с помощью моделирования тестов (программа Т01ЮН2) путем подбора проницаемости во всей расчетной области. В результате получена трехмерная карта проницаемости на сетке 34x26x24.

Эти данные не могут, конечно, характеризовать структуру сети трещин на исследуемом участке. В то же время они дают информацию о разбросе локальных значений проницаемости. Этот разброс очень велик — пять порядков величины. Быстрый отклик давления, наблюдаемый во всех скважинах, свидетельствует о том, что все они хорошо соединены друг с другом сетью трещин (по крайней мере, пневматически).

Многочисленные исследования показывают, что стандартные гомогенизированные модели приводят к качественно неверным результатам. Причиной таких отклонений может служить фрактальный характер сетей трещин, по которым происходит просачивание влаги, изменяющий сам характер осредненного закона Дарси.

Эту задачу можно формулировать как двухфазную, когда первой фазой является вода, а второй - вытесняемый ей воздух. Такая постановка приводит к системе из двух уравнений Лаверетта-Бакли, описывающих баланс каждой фазы. В случае, когда второй фазой является нефть, система уравнений Лаверетта-Бакли достаточно хорошо изучена и широко применяется на практике. Плотности и подвижности воды и нефти относительно близки, поэтому трудностей принципиального характера при численном решении этих уравнений не возникает, либо они достаточно легко преодолеваются. Иное дело, когда второй фазой является воздух или газ, отличающийся по этим параметрам от воды на несколько порядков.

Расчеты задач просачивания в сильно неоднородных трещиновато-пористых средах представляют собой серьезную вычислительную проблему, поскольку наличие у большинства известных алгоритмов аппроксимационной вязкости (в рассматриваемом случае «аппроксимационных капиллярных сил») существенно искажает характер решения в экстремальных случаях. Естественным решением, в этой ситуации, представляется выбор численной схемы, относящейся к т.н. алгоритмам «высокой разрешающей способности» [15, 17, 18, 21, 24, 25, 33]. Альтернативным подходом можно считать использование схемы КАБАРЕ [34-45]

Схема КАБАРЕ для простейшего одномерного линейного уравнения переноса была предложена и подробно исследована в работах Головизнина В.М. и Самарского A.A. в 1998 году. Позже выяснилось, что в западной литературе она известна как схема Айзерлиса (Upwind Leapfrog) поскольку является представителем семейства разностных схем, исследованных эти автором на устойчивость в 1986 году. В дальнейшем, эта схема претерпела ряд эволюционных скачков. Значительный вклад в ее развитие внесли Головизнин В.М., Карабасов С.А., Кобринский И.М. К наиболее важным изменениям исходной трехслойной схемы КАБАРЕ (Upwind Leapfrog) можно отнести ее представление в двухслойном виде, что было осуществлено введением дополнительных неизвестных, т.н. «консервативных» переменных. Затем она была дополнена алгоритмом коррекции потоков, базирующемся на прямом использовании принципа максимума. Это сделало возможным ее обобщение на более содержательные и сложные по сравнению с простейшим уравнением переноса, случаи.

Для задач подземной гидродинамики схема КАБАРЕ была впервые использована в работе (Головизнин В.М., Карабасов С.А.,1998) для решения системы уравнений Лаверетта-Бакли в двумерной задаче о скважине. Работа носила чисто методический характер - в ней было проведено сравнение четырех возможных способов обобщения схемы КАБАРЕ, три из которых основывались на трехслойной версии алгоритма. Результат, полученный по одной из этих схем, оказался на то время удовлетворительным. Дальнейшее развитие схемы КАБАРЕ привело к необходимости критического пересмотра полученных ранее результатов с целью их качественного улучшения на новой методической основе.

Целью настоящей работы является^дальнейшее развитие схемы КАБАРЕ для решения одномерных нелинейных законов сохранения с выпуклыми и невыпуклыми функциями потоков; разработка на базе модифицированной схемы нового эффективного вычислительного алгоритма для численного решения одномерных задач однофазной и двухфазной фильтрации в зонах неполного влагонасыщения; разработка новых математических моделей двумерной однофазной и двухфазной фильтрации в геологических формациях с сильно неоднородными фильтрационными свойствами.

Диссертация состоит из введения, трех глав основного текста, заключения и списка литературы.

Заключение

Разработан новый вид численных потоков для схемы КАБАРЕ. Проведены исследования пределов применимости новой модификации методики КАБАРЕ для задач с невыпуклой функцией потоков. Показано, что включение в алгоритм диссипатора Паниковского с коэффициентом £ = 0.6 позволяет получать энтропийные решения для невыпуклых функций потока при высоких порядках касания производной в точках смены выпуклости. Установлено, что к такому же эффекту приводит сужение области возможных значений потоковых переменных, допускаемой принципом максимума.

На основе схемы КАБАРЕ разработаны новые вычислительные алгоритмы для одномерной однофазной и двухфазной фильтрации в зоне неполного влагонасыщения. Задача однофазной фильтрации сформулирована как задача решения одномерного квазилинейного гиперболического уравнения при наличии ограничивающей связи.

Разработан алгоритм учета такой связи в алгоритме КАБАРЕ с использованием вариационного принципа для возмущений потоков. Предложен способ учета капиллярных эффектов. Для двухфазной фильтрации с учетом явления вытеснения воздуха получено аналитическое выражение для градиента давления через величину влагонасыщения, что привело к классической одномерной гиперболической задаче с немонотонной функцией потоков. Проведено сравнение численных результатов решения модельной задачи о просачивании по двум этим моделям, Показано, что при реальных значениях отношения подвижностей воздуха и воды физические результаты получаются близкими, в то время как вычислительные затраты в однофазной модели существенно меньше.

Разработаны новые математические модели двумерной однофазной и двухфазной фильтрации в зоне аэрации. Для двумерной двухфазной модели Лаверетта-Бакли с учетом вытесняемого воздуха разработан вычислительный алгоритм на основе схемы КАБАРЕ. Вычислительные особенности этого алгоритма продемонстрированы на модельной задаче о напорном течении в пористой среде с непроницаемой трещиной.

Представлена двумерная постановка задачи об однофазной фильтрации с учетом ограничений на решение. Разработан эффективный вычислительный алгоритм, являющийся обобщением одномерного алгоритма. Для убыстрения счета задач о динамике протечек на перколяционных решетках разработано т.н. «гидростатическое» приближение. Работоспособность и вычислительная эффективность нового подхода продемонстрирована на ряде модельных и тестовых задач.

1. Gwo, J.P., et al., Using multiregion model to study the effect of advective and diffusive mass transfer on local physical nonequilibrium and solute mobility in a structured soil. Water Resources Research, 1996. 32(3): p. 561-570.

2. Pruess, K., A mechanistic model water seepage through thick unsaturated zones in fractured rocks of low matrix permeability. Water Resources Research, 1999. 35(4): p. 1039-1051.

3. Баренблатт, Г.И., B.M. Ентов, and B.M. Рыжик, Движение жидкостей и газов в пористых пластах. Недра, 1984: р. 224.

4. Коновалов, А.Н., Задачи фильтрации многофазной несжимаемой жидкости. «НАУКА», Сибирское отделение, 1988: р. 166.

5. Самсонов, Б.Г. and JI.M. Самсонова, Миграция вещества и решение гидро-геологических задач. Недра, 1987: р. 118.

6. Селяков, В.И. and В.В. Кадет, Перколяционные модели процессов переноса в микронеоднородных средах. Недра, 1995: р. 222.

7. Закиров, С.Н., et al., Многомерная и многокомпонентная фильтрация. 1988, Москва: Недра. 335.

8. Коллинз, Ю.Р., Течение жидкостей через пористые материалы, ed. Г.И. Баренблат. 1954, Москва.

9. Нигматулин, ¥ Ж., Динамика многофазных сред. 1987, Москва: Наука. 464,360.

10. Николаевский, В.Н., е1 а1., Механика насыщенных пористых сред. 1970, Москва: Недра. 339.

11. Цыпкин, Г.Г., Аналитическое решение нелинейной задачи разложения газового гидрата в пласте. Механика жидкости и газа, 2007. 5: р. 133142.

12. Цыпкин, Г.Г., Инжекция раствора соли в геотермальный резервуар, насыщенный перегретым паром. Механика жидкости и газа, 2008. 5: р. 120-131.

13. Цыпкин, Г.Г., Влияние капиллярных сил на распределениевлагонасыщенности при протаивании мерзлого грунта. Механика жидкости и газа, 2010. 6: р. 122-132.

14. Цыпкин, Г.Г., Нелинейная задача протаивания ненасыщенного мерзлого грунта при наличии капиллярных сил. Механика жидкости и газа, 2012. 1.

15. Boris, J.P. and D.L. Book, Flux-corrected transport I. SHASTA, a fluid transport algorithm that works. J. Comput. Phys., 1997. 135(2): p. 172-186.

16. Fedkiw, R.P., B. Merriman, and S. Osher, Efficient characteristic projection in upwind difference schemes for hyperbolic systems: the complementary projection method. J. Comput. Phys., 1998. 141(1): p. 22-36.

17. Harten, A., et al., Uniformly high order accurate essentially non-oscillatory schemes, 111. J. Comput. Phys., 1987. 71(2): p. 231-303.

18. Jiang, G.-S., et al., High-Resolution Nonoscillatory Central Schemes with Nonstaggered Grids for Hyperbolic Conservation Laws. SIAM J. Numer. Anal., 1998. 35(6): p. 2147-2168.

19. Kurganov, A. and E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection ¡diffusion equations. J. Comput. Phys., 2000. 160(1): p. 241-282.

20. Li, Y., Wavenumber-extended high-order upwind-biased finite-difference schemes for convective scalar transport. J. Comput. Phys., 1997. 133(2): p. 235-255.

21. Marano, S. and M. Franceschetti, Ray propagation in a random lattice: a maximum entropy, anomalous diffusion Process. IEEE Transactions on Antennas and Propagation, 2005. 53(6): p. 1888-1896.

22. Mazzia, A., L. Bergamaschi, and M. Putti, A time-splitting technique for the advection-dispersion equation in groundwater. J. Comput. Phys., 2000. 157(1): p. 181-198.

23. Mohanty, B.P., et al., New piecewise-continuous hydraulic functions for modeling preferential flow in intermitten-flood-arranged field. Water Resources Research, 1997. 33(9): p. 15.

24. Sheu, T.W.H., S.K. Wang, and S.F. Tsai, Development of a high-resolution scheme for a multi-dimensional advection-diffusion equation. J. Comput. Phys., 1998. 144(1): p. 1-16.

25. Shu, C.-W. and S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes,II. J. Comput. Phys., 1989. 83(1): p. 3278.

26. Yavneh, I., Analysis of a fourth-order compact scheme for convection-diffusion. J. Comput. Phys., 1997.133(2): p. 361-364.

27. Peters, R.R. and E.A. Klavetter, A continuum model for water movement in an unsaturated fractured rock mass. Water Resour. Res., 1988. 24(3): p. 416430.

28. Wang, J.S.Y. and T.N. Narasimhan, Hydrologic Mechanisms Governing Fluid Flow in a Partially Saturated, Fractured, Porous Medium. Water Resour. Res., 1985. 21(12): p. 1861-1874.

29. Pruess, К., B. Faybishenko, and G.S. Bodvarsson, Alternative concepts and approaches for modeling flow and transport in thick unsaturated zones of fractured rocks. Journal of Contaminant Hydrology, 1999. 38(1-3): p. 281322.

30. Dahan, O., et al., Field observation offlow in a fracture intersecting unsaturated chalk. Water Resour. Res., 1999. 35(11): p. 3315-3326.

31. Dahan, O., et al., On Fracture Structure and Preferential Flow in Unsaturated Chalk. Ground Water, 2000. 38(3): p. 444-451.

32. Huang, K., Y.W. Tsang, and G.S. Bodvarsson, Simultaneous inversion of air-injection tests in fractured unsaturated tuff at Yucca Mountain. Water Resour. Res., 1999. 35(8): p. 2375-2386.

33. Куликовский, А.Г., H.B. Погорелов, and А.Ю. Семенов, Математические вопросы численного решения гиперболических систем уравнений. 2001, Москва: Физматлит. 607.

34. Goloviznin, V.M., A balance-characteristic method for the numerical solution of one-dimensional equations of gas dynamics in Euler variables. Mat. Model., 2006.18(11): p. 14-30.

35. Goloviznin, V.M. and S.A. Karabasov, Digital Transport Algorithm for Hyperbolic Equations. Hyperbolic problems: theory, numerics, applications. 2006, Yokohama: Yokohama Publishers.

36. Goloviznin, V.M. and S.A. Karabasov, Compact Accurately Boundary-Adjusting High-Resolution Technique for Fluid Dynamics. J. Comput. Phys., 2009. 228(19): p. 7426-7451.

37. Goloviznin, V.M., S.A. Karabasov, and I.M. Kobrinskii, Balance-characteristic schemes with separated conservative and flux variables. Mat. Model., 2003.15(9): p. 29-48.

38. Goloviznin, V.M. and A.A. Samarskii, Some properties of the difference scheme "cabaret". Mat. Model., 1998. 10(1): p. 101-116.

39. Goloviznin, V.M. and A.A. Samarskii, Finite approximation of convective transport with a space splitting time derivative. Mat. Model., 1998. 10(1): p. 86-100.

40. Головизнин, B.M., Балансно-характеристический метод численного решения уравнений газовой динамики. Докл.Акад. Наук, 2005. 403(4): р. 1-6.

41. Головизнин, В.М. and С.А. Карабасов, Дискретные математические модели двухфазной фильтрации с пространственным расщеплением временной производной. Известия РАН, серия Энергетика, 2000(4).

42. Головизнин, В.М. and С.А. Карабасов, Некоторые примеры численного моделирования двумерной фильтрации. Препринт ИБРАЭ. 1998, Москва: ИБРАЭ РАН.

43. Головизнин, В.М. and С.А. Карабасов, Нелинейная коррекция схемы «КАБАРЕ» Математическое Моделирование, 1998. 12(1): р. 107-123.

44. Головизнин, В.М., et al., Новый вычислительный алгоритм для математического моделирования просачивания влаги сквозь ненасыщенную трещиноватую геологическую среду с низкой проницаемостью. Препринт ИБРАЭ №IBRAE 2006-07. 2006, Москва: ИБРАЭ РАН. 53.

45. Iserles, A., Generalized Leapfrog Methods. IMA Journal of Numerical Analysis, 1986. 6: p. 381-392.

46. Ostapenko, V.V., On the monotonicity of the balance-characteristic scheme. Mat. Model., 2009. 21(7): p. 29-42.

47. Горицкий, Ф.Ю., C.H. Кружков, and Г.А. Чечкин, Уравнения с частными производными первого порядка. 1999, Москва: МГУ им. М.В. Ломоносова (учебное пособие). 95 стр.

48. Kurganov, A., G. Petrova, and В. Popov, Adaptive Semidiscrete Central-Upwind Schemes for Nonconvex Hyperbolic Conservation Laws. SIAM J. Sci. Comput., 2007. 29(6): p. 2381-2401.

49. Kurganov, A., S. Noelle, and G. Petrova, Semidiscrete Central-Upwind Schemes for Hyperbolic Conservation Laws and Hamilton—Jacobi Equations. SIAM J. Sci. Comput., 2001. 23(3): p. 707-740.

50. Pruess, К., C.M. Oldenburg, and G.J. Moridis, TOUGH2 User's Guide Version 2. 1999, Berkeley: Lawrence Berkeley National Laboratory.

51. Buckley, S.E. and M.C. Leverett, Mechanism of Fluid Displacement in Sands. Petroleum Transactions, AIME, 1942.146: p. 10.

52. Fischer, U., et al., Modeling nonwetting-phase relative permeability accounting for a discontinuous nonwetting phase. Soil Science Society of America, 1997. 61(5): p. 15.

53. Ортега, Д. and P. Вернер, Итерационные методы решения нелинейных систем уравнений со многими неизвестными. 1975, Москва: МИР.

54. Finsterle, S., J.Т. Fabryka-Martin, and J.S.Y. Wang, Migration of a water pulse through fractured porous media. Journal of Contaminant Hydrology, 2002. 54(1-2): p. 37-57.

55. Harten, A., High resolution schemes for hyperbolic conservation laws. J. Comput. Phys., 1997. 135(2): p. 260-278.

56. Tsypkin, G.G., Mathematical Models of Gas Hydrates Dissociation in Porous Media. Annals of the New York Academy of Sciences, 2000. 912(1): p. 428436.

57. Басниев, K.C., И.Н. Кочина, and B.M. Максимов, Подземная гидромеханика. 1993: Недра. 417.

58. Лукхаус, С. and П.И. Плотников, Энтропийные решения уравнений Баклея-Леверетта. Сибирский математический журнал, 2000. 41(2): р. 400-420.