Метод смешанной переменной для моделирования насыщенно-ненасыщенных течений Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

Том 149, ки. 4

Физико-математические пауки

2007

УДК 517.9

МЕТОД СМЕШАННОЙ ПЕРЕМЕННОЙ ДЛЯ МОДЕЛИРОВАНИЯ НАСЫЩЕННО-НЕНАСЫЩЕННЫХ ТЕЧЕНИЙ

A.A. Ахтареев, Р.З. Даутов

Аннотация

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

Введение

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

дв

ф— = х € О, (1)

ве = Б(ф). (2)

Здесь д - поток жидкости, д = —К(ве)К0(Уф + в), И С Яй - область течения, занятая пористой средой, й = 1, 2, 3; эффективная насыщенность в е связана с в соотношением

в в*

Параметры пористой среды в* и в*, 0 < в* < в* < 1, называются соответственно остаточной и предельной насыщенностями; ф - пористость среды (ф £ (0,1)); Ко = Ко(х) — тензор проводимости насыщенной среды. Постоянный единичный вектор в равен орту оси, противоположной направлению силы тяжести (в = 0, если силой тяжести можно пренебречь), функция Q описывает источники и стоки жидкости. Пористая среда называется насыщенной, если ве = 1; ненасыщенной, если 0 < ве < 1; сухой, ее ли ве = 0.

Функции К (относительная проводимость среды) и Б также являются параметрами пористой среды и считаются заданными. Функция Б(ф) монотонно возрастает на (—те, го), Б(ф) = 1 при ф > фа, Б(ф) ^ ^и ф ^ —го; К(ве) является

Рис. 1. Аппроксимации Ван-Генухтена функций K(se) и щи а = 1

монотонно возрастающей и выпуклой, К(0) = 0, К(1) = 1. Наиболее распространены аппроксимации Брукса Кори [1] и Ван-Генухтена [2] этих функций. Аппроксимации Брукса Кори (ВС) имеют вид

КЫ = адЛ1' ф ~-1/а

К ' е ^ ' \(—аф)-", ф< — 1/а

(то есть фа = —1/а). Согласно Ван-Генухтену (УС) имеем:

K(se) = s1/2 (l - (1 - 2 , S(ф)=|

1, ф > 0,

(1 + (—аф)п)-т , ф< 0.

Здесь п > 1, а > 0 - параметры среды, т = 1 — 1/п. На рис. 1 приведены графики этих функций для аппроксимации УС при двух значениях п. Отметим, что в этом случае К(ве) ведет себя как в0'5+2/т при малых ве; Б(ф) ~ 1/(—ф)п-1 при ф —то; Б близка к ступенчатой функции при больших п.

Пористая среда называется однородной, если все параметры среды не зависят от х облает и О, в противном случае - неоднородной. Далее мы ограничимся рассмотрением только кусочно-однородных пористых сред. А именно, будем предполагать заданным разбиение области О та совокупность подобластей О^ таких, что параметры среды в них постоянны, кроме, возможно, функции Ко. Для таких сред к уравнениям (1), (2) присоединяются соотношения непрерывности функции ф и нормальной составляющей потока д та общих границах областей О^. Отметим, что в случае неоднородной среды в уравнениях (1), (2) необходимо полагать К(х'ве), Б(х, ф) вместо К(ве) и Б(ф) соответственно. Однако для краткости мы

х

Уравнения (1), (2) дополним следующими краевыми и начальным условиями:

ф = фо, х е Го, — д ■ V = до, х е гь г > 0' (3)

ф = ф0, х е О, г = 0. (4)

Здесь Г = Го и Г1 — граница О, V — единичный вектор внешней нормали к Г.

Уравнение (1) с учетом (2) представляет собой нелинейное уравнение в частных производных для определения давления. Оно является вырождающимся и имеет переменный тип. Действительно, в тех точках области течения, где ве = 0 (в = = в*), то есть в сухой области, оператор ф ^ д является вырождающимся, и,

следовательно, вырождается левая часть уравнения (1). В области {0 < ве < 1} неполного насыщения уравнение (1) относится к параболическому типу, в области {ве = 1} полного насыщения - к эллиптическому.

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

1. Известные подходы к построению численных методов решения задачи

Кратко опишем широко распространенные в практике вычислений подходы к

аппроксимации уравнений (1), (2), отмечая их основные достоинства и недостатки. ф

ную насыщенности, придем к одному уравнению

Ф'Щ^р- - <ПчЦФ)Ко(УФ + е) = д, X € п, (5)

где S('ф) = Б (ф), К(ф) = К (Б (ф)), фя = ф(в* — в*). Далее, на основе метода конечных разностей или метода конечных элементов это уравнение аппроксимируется и решается каким-либо итерационным методом (обычно методом Ныотона или его вариантами). По найденному ф неизвестная в находится из уравнения (2). Для аппроксимации по времени наиболее часто используется неявный метод Эйлера и разностные схемы второго порядка точности, которые могут комбинироваться с соответствующими процедурами автоматического выбора шага интегрирования.

К несомненным достоинствам этого подхода, кроме сокращения числа неизвестных вдвое, следует отнести его универсальность: метод применим к насыщенно-

ненасыщенным неоднородным средам произвольной структуры. Основной недо-

ф

подобластей малой насыщенности (достаточно сухих зон). Поясним сказанное подробнее.

При наличии в среде достаточно сухих зон, где ве = 0 или ве « 0, решение задачи имеет ярко выраженный фронтовой характер, причем фронт тем круче, чем ближе функция Б к ступенчатой. При этом на фронте ф « —го, |Уф| « го, К(ф)Уф « 0. Таким образом, неизвестная ф плохо определена в окрестности фрон-—го

свойство в большей степени является отрицательным для итерационных методов, требующих хорошего начального приближения, поскольку в этом случае решение в момент времени £ является плохим приближением к решению при £ + Д£. Например, при использовании итерационных методов типа Ныотона это приводит к достаточно сильному ограничению на шаг сетки по времени Кроме того, наличие сухих зои приводит к вырождению левой части уравнения (5) в этих зонах, что существенно ухудшает обусловленность сеточных уравнений и замедляет сходимость итерационных методов. Напротив, наличие в области течения зон полного

насыщения не вызывает существенных трудностей при данном подходе.

ф

Ф

-&у{К0ЧО,{ф)+1С{ф)К0е) = д, 0.(Ф) = I К{Б{ф))с1ф, (6)

о

аппроксимация которой приводит к более простым по структуре сеточным уравнениям. Отметим также, что эта форма позволяет ввести новую неизвестную u = - потенциал Кирхгофа, в терминах которой задача приобретает вид

- div (A'oVw, + К(ы)ВД = Q,

где S(u) = S(Q-1(u)), K(u) = K(Q-1(u)), S'(u) ^ щи u ^ +0. Новая неизвестная u является лучше определенной, чем ф, а уравнение более не является вырождающимся по пространственным переменным. Трудности возникают в тех

S

зования Кирхгофа является то, что оно применимо только к однородным средам (в отличие от уравнений (5), (6)).

Методы на основе s -формы уравнений. Если исходные данные таковы, что зоны полного насыщения не реализуются, то есть se(x, i) < 1 при любых x, i, то в уравнениях (1), (2) можно исключить неизвестную ф. Придем к уравнению

ds

<t>a-:£-àiv(KoD(se)Vse + K(se)Koe) = Q, D(se) = K(se)P'(se), dt

или

Se

ds Г

-div(K0VQ(se)+K(se)K0e) =Q, Q(se) = D{s)ds,

(7)

0

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

при se ^ 0 : D(se

|2+1/n, для ВС, I 0.5 + m, для VG.

Однако возникают трудности при наличии в области течения зон, близких к полному насыщению, что проявляется в поведении функции диффузивности П: П(ве) ^ при ве ^ 1. Основной недостаток данного подхода заключается в том, что он применим только к однородным ненасыщенным средам (для неоднородных сред функция ве является разрывной на общих границах сред). При решении уравнения (7) также можно ввести преобразование Кирхгофа.

Методы на основе смешанной формы уравнений. Методы этого типа основаны на непосредственной аппроксимации уравнений (1), (2) с дальнейшим исключением одной из неизвестных на сеточном уровне на итерационном шаге. К новым схемам решения приходим в том случае, если такое исключение осуществлять динамически (в зависимости от состояния среды). Опишем один такой метод, ограничиваясь случаем однородной среды [4].

Уравнение (1) аппроксимируем неявным методом конечных разностей или конечных элементов. В каждом узле сетки ж;, I = 1,..., и, определены две неизвестные в; = ве(ж;) и ф; = ф(ж;), но только одна из них независимая. Какую неизвестную считать первичной, а какую вторичной (то есть зависимой), определяется нашим желанием и обратимостью соотношения (2). Сеточную аппроксимацию в момент времени £ + условно запишем в виде А;(в, ф) = 0, I = 1,..., и.

Обозначим первичную переменную в узле ж; через и;. Тогда сеточные уравнения примут вид А; (и) = А; (в, ф) = 0, I = 1,...,и. Эти уравнения решаются при помощи метода Ньютона 7(ик)Дик+1 = — А(ик), Дик+1 = ик+1 — ик, при

следующим определении первичных переменных на каждой итерации:

если вгк < ва, если вгк > вь, иначе.

Последний случай означает, что принятый способ определения первичной переменной сохраняется таким же, как и на итерации. 3десь 0 < ва < вь < 1 - параметры метода (например, ва = 0.89, вь = 0.99). Способ вычисления матрицы Якоби ,/(ик) следующий: /-я строка J(ик) равна /-й строке матрицы дА(вк,фк)/дв, если игк+1 есть в переменная, иначе равна /-той строке матрицы дА(вк,фк)/дф.

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

к+1

к+1

Т1, к,

2. Новая неизвестная. Эквивалентная переформулировка задачи для однородной среды

Приведем новую формулировку исходной задачи (1), (2), которую в случае однородных сред можно рассматривать как непрерывный аналог описанного выше метода первичных переменных.

Ранее отмечалось, что в зоне малых насыщенностей, например, при ф < фа,

хорошо определена неизвестная ве, а в зоне больших насыщенностей - неизвестная ф

и = и(ф)={^(Ф)' Ф < ' (8)

[Б'(фа)(ф - фст) + Б(фа), ф>фст,

обратная к которой равна

ф = Р(и)={Р(и)' , " < ^)' (9)

\(и - Б(фа))/Б'(фст) + фа, и>Б(фа).

По определению функция и = и(ж,£) является «смешанной» (составной) неизвестной: если и(ж,£) < ва, ва = Б(фа), то и(ж,£) совпадает с эффективной насыщенностью ве(ж,£), в противном случае — с «нормированным» давлением Б'(фа)(ф — Фа) + ва. По известной и давление вычисляется по формуле (9), а насыщенность определяется по формуле ве = Б(ф).

фа

гиба функции Б(ф). В этом случае функции и = и(ф) и ф = Р(и) являются дважды непрерывно дифференцируемыми (иначе непрерывно-дифференцируемыми). Нетрудно подсчитать, что ва = (1/(1 + ш))т. Для аппроксимации ВС можно взять ва € [0.8,0.9].

Введем следующие функции переменной и € [0, го):

5(и) = I и' и < ва,

(Б ((и — ва )/Б'(фа )+ фа) ' и > ва,

и

К(и) = К(5(и)), Р(и) = К(5(и))Р'(и), 2(и) = Р(и) ¿и

0

Тогда исходные уравнения (1). (2) перепишутся в виде системы

- <Цу {К0К(зе)(ЧГ(и) + е)) = д,

- <Ну (К0Ци)(УГ(и) +е))= д.

Последнее удобно записать в следующем виде:

(Ю)

Это уравнение дополняется краевыми и начальным условиями (см. (3). (4)):

Здесь 9 = —Ко(У2(и) + К(и)е).

Задача (10). (11) представляет собой искомую переформулировку исходной задачи в терминах смешанной переменной. Отмстим, что в тех точках пористой среды, где и < ва, уравнение (10) совпадает с в-формой уравнения Ричардса, иначе

ф

ства обеих формулировок. После решения задачи (10), (11) определяем решение исходной задачи по формулам

Отметим также формальное тождество уравнений (6) и (10). Поэтому для их аппроксимации и решения полученных сеточных уравнений могут использоваться одни и те же численные методы. Разница между этими уравнениями заключается в свойствах их коэффициентов и самой неизвестной, что существенно влияет на эффективность приближенных методов. Для сравнения графики коэффициентов 5, <2, К уравнений (6) и (10) для аппроксимации Ув при и =10 представлены на рис. 2. Па оси абсцисс отмечены значения фст, ва и иа = ва — Б'(фст)фст,

ф=0

Анализ приведенных функциональных зависимостей, а также графиков, позволяет сделать следующие выводы о различии уравнений (6) и (10) в случае ненасыщенной фильтрации.

1. Переменная ф меняется в неограниченном интервале (—го, фа], и - в ограниченном (0, иа].

2. Уравнение (6) вырождается по всем коэффициентам в сухой зоне, уравнение (10) нет. Главное различие коэффициентов уравнений - в свойствах функций 5(ф)

и 5(и): 5(и) линейна при малых и, тогда как Б(ф) ^ 0 при больших отрицатель-ф

3. Функция Б(ф), в отличие от 5(и), может быть близка к ступенчатой функции при определенных параметрах среды.

Отмеченные особенности уравнения (10), сочетающего в себе достоинства в- и

ф

шения и, соответственно, для моделирования насыщенно-ненасыщенных течений в пористых средах.

и = ио = и(фо), ж € Го, —9 • V = 90, ж € Г1, г > 0, и = и0 = и(ф0), ж € О, г = 0.

(Н)

ф = Р (и),

5(и).

Рис. 2. Коэффициенты S, С, Q уравнения VG при n = 10 а = Ко = фв = 1

u

(слева) и (10) (справа) для аппроксимации

3. Случай неоднородной среды. Вычисление коэффициентов уравнения

Описанная выше методика замены неизвестных непосредственно обобщается на случай кусочно неоднородных пористых сред.

Пусть известно разбиение области О та совокупность подобластей (сред) Oj таких, что параметры среды в них постоянны, кроме, возможно, функции Ко. Пусть щ, aj - параметры VG или ВС для г-й среды О^ se = Sj(^), ф = Pj(se) -капиллярные зависимости, Kj(se) - относительная проводимость г-й среды, определяемые по nj и aj согласно аппроксимациям VG или ВС. Обозначим

n = nk = min nj, а = a^. j

Неизвестную u определим то формулам (8), (9) для параметров k-й среды, положим

ф = P() =! Pk(u), u < ska,

kW \(u - 4)^(ф£)+ фка, u>sk . Определим также функции

Sj(u) = Sj (Pk (u)), Kj(u) = Kj(Sk(u)),

u

Dj(u) = Kj(u)Pk(u), Qj(u) = У Dj(u) du.

о

В новых переменных задача по-прежнему примет вид уравнений (10), (11), в которых теперь функции S, Q, K равны S^ Qj; ^ в г-й области О^ Ф«|п- = = ^j(s* — s*,j). Несложный анализ функции S(u) показывает, что уравнение (10) остается вырождающимся в тех подобластях О^, для которых nj = n. Например, если S(ф) ~ (—ф ^ —го в О^ то S(u) ~ uni/n при малых u в Oj. Таким

n

раз меньше.

Функции S, K, Q достаточно сложно определяются и требуют немалого числа

Q

Сделаем практическое замечание, касающееся способа вычисления как самих этих

функций, так и их производных. Учитывая, что нетривиальное изменение функций имеет место лишь на небольшом отрезке [0, иа], иа = 0(1), то естественно перейти на этом отрезке от функций 5, К, 2 к их представлениям сплайнами. Далее, мы используем непрерывную кусочно-линейную интерполяцию на равномерной сетке с шагом Ни = иа/(пи — 1). Для полноты изложения приведем следующий хорошо известный алгоритм.

1. В узлах щ = (г — 1)^и вычисляем значения 5(и»), К(и»), Р(и»), г = 1,..., пи, для каждой зоны . Для каждой зоны также вычислявм значения 2 (и») = = 2(и»_1) + 0.5(Р(иг_1) + Р(и»))Л.и, г = 2,..., пи, 2(и^ = 0. Это приближения интеграла по формуле трапеций.

Пусть теперь / означат любую из функций 5, ^ шш 2.

2. Вычисляем с/(и») = (/(и»+1) — /(и»))/^, г = 1,..., Пи — 1.

Все эти вычисления проводятся один раз до начала решения задачи.

3. Для вычисления /(и) и /'(и) используем следующие операции:

- найдем отрезок [и», и»+1], содержащий точку и: г = /1оог(и/^и + 1);

— если и < иа положим /(и) = /(и») + с/(и»)(и — и»), /'(и) = с/(и»);

Вне интервала [0, иа] функции 5 и К постоянны и равны единице, а функция 2 линейна. Этот метод за минимум арифметических операций обеспечивает точность вычисления функций 5, 2, К порядка 10_6 ^ 10_8 щи пи порядка 103 ^ 104. Этого вполне достаточно для практических целей, если учесть, что параметры среды известны также приближенно.

4. Аппроксимация задачи в новых переменных

Пусть [0, Т] — интервад изменения переменной Введем на нем неравномерную сетку узлов 0 = ¿о < ¿1 < ... < ¿м = Т с шагами т^ = ^ —• Для аппроксимации по времени задачи (10), (11) используем неявную схему Эйлера

- «ну {Коча(и) + ¡с(и)Кое) = д,

и = ио, х £ Го, —9 • V = ^о, х € Г1.

Здесь все функции вычисляются в текущий момент времени Ь = Ь^, кроме и = = и(х, ¿¿_1), т = т^-, = 1,..., М. В начальный момент времени решение известно: и| о = и0. Для аппроксимации этой эллиптической краевой задачи используем метод конечных элементов на основе линейных треугольных или билинейных четырехугольных изоиараметрических конечных элементов. В дополнение к стандартному методу используем две дополнительные аппроксимации, которые не понижают теоретической точности метода.

1. Аппроксимируем функции 5, 2, К в слабой формулировке задачи так же, как и и. Точнее, пусть щ — значение функции и в узле сетки х», уч(х) — базисная функция, соответствующая этому узлу, г = 1,..., N. Тогда аппроксимации и и 5 определяются как

N N

»=1 »=1

Аналогично определяются К^(и^) и 2^(и^).

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

Тогда

/л

f(x) dx = J2 J f(x)dx « E — E le l/(-<) =

П e e e e i=1

Здесь |e| - площадь e; же, i = 1,..., me - координаты вершин элемента e (me = 3 или me = 4). Эта квадратура позволяет диагонализировать матрицу масс. Аналогично для вычисления интеграла по границе Г1 используем составную формулу трапеций Sri.

Искомая дискретная задача определяется тождеством:

- 5л(йл))/г + (KoVQh(uh) + Kh(uh)Koe) • Vvh) =

= Sn(Qvh) + Sri (qovh),

где пробная функция vh предполагается равной нулю на Г0.

Использование указанных аппроксимаций функций S, Q, K позволяет существенно упростить как вычисление невязки сеточного уравнения (R), так и его матрицы Якоби (J). Например, в случае однородной среды имеем (посредством элементной сборки):

R(u) = D(S(u) - S(u)) + AQ(u) + CK(u) - F = 0, (12)

J (u) = DS'(u) + AQ'(u) + CK'(u).

Здесь u = (u1;..., uN)T - вектор узловых параметров; f (u) = (f (u1),..., f (uN))T, f'(u) = diag(f '(ui),..., f'(uN)) - диагональная матрица, f = S, Q, K; диагональная матрица D и матрицы A C являются конечно-элементными аппроксимациями (с численным интегрированием) соответственно линейных операторов

—div (Ä'nV(-)), -div (Ä'ne(-)).

В случае неоднородной среды считаем, что разбиение на элементы согласовано с границами подобластей, а конечный элемент полностью относится к одной подобласти. Используется также элементная сборка R и J. Укажем, например, способ вычисления матрицы AQ'(u) (слагаемо ев J (u)). Остальные составляющие R(u) J(u) e

ui, uf,..., um,e - его узловые параметры (то есть компоненты u, относящиеся к вершинам e), Ae - м атрица жесткости эле мента e, соответствующая линейному оператору -div (K0V(-)), Q'(ue) = diag(Q'(ui),..., Q'(ume))• Тогда AQ'(u) получается стандартной процедурой суммирования локальных матриц AeQ'(ue). Отметим, что это совсем дешевая процедура на итерациях, если хранить постоянные матрицы Ae для всех конечных элементов.

5. Решение сеточных уравнений.

Управление шагом интегрирования

Хорошо известно, что эффективность метода Ныотона при решении нестационарных сеточных задач критическим образом зависит от выбора шагов интегрирования по времени. Если шаг интегрирования т велик, то метод Ньютона либо не сходится, либо сходится медленно. Напротив, если шаг т мал, то достаточно выполнить лишь одну итерацию. Известны различные стратегии автоматического выбора шага, целыо которых является достижение максимальной эффективности и точности всего процесса интегрирования. При решении тестовых задач мы использовали алгоритм из [3], хорошо зарекомендовавший себя при моделировании

течений в пористых средах [4. 5]. Этот метод использует единственный параметр 6 - допустимую ошибку дискретизации по времени как для контроля сходимости

т

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

В каждый момент времени вычисляется решение и и вектор и.

0. Начальный шаг. При Ь = 0 вектор и находится как решение при и = ио сеточного аналога задачи

ф85'(и)и = (КоУ2(и) + К(и)Кое) + д, и = 0, х € Го, д • V = ^о, х € Г1.

т1

Пусть при Ь = решение и, и и т = т, известны и задан параметр погрешности 6. Для определения решения и при Ь = Ь, выполняются следующие шаги:

1. Шаг (предиктор). Определяем ир = и + та;

2. Шаг (корректор). Делаем один шаг по методу Ныотона с начальным приближением ир и находим и, норму поправки С = 0.25||и — ир||2/||и||то, а также и = (и — и)/т. Здесь

1 N 1/2

1М1 2 = ) ' 1Н1 = ^^

' =1

3. Выбор шага интегрирования. Полагаем тр = т(6/С)°-5. Если тр > т, то шаг интегрирования считается выбранным успешно и полагается т,+1 = = тш{ттах,дт, тр}, д = 2,3. Если (г < тр < т (С € [0.8,0.9]), то решение и

т,+ 1 = т. В противном случае, шаг интегрирования считается выбранным неудачно и вычисления повторяются с сокращенным шагом т = тах{тт^,т26/(Стр)}.

Для решения несимметричной разреженной системы линейных алгебраических уравнений (12) на итерации Ныотона

7(ир)(и — ир) = —Д(ир)

нами использовался итерационный метод ВЮС^аЬ [6] с неполным модифицированным Ьи-разложением матрицы в качестве предобуславливателя (МПДДО)) [7]. Выход из итераций осуществлялся при выполнении условия

норма невязки < 10_6

норма правой части

е=0

задаче с симметричной матрицей после введения новой неизвестной Д = 2'(ир)(и— —ир). В этом случае вместо BiCGstab достаточно использовать метод сопряженных градиентов.

Указанный выбор стратегии управления шагом и метода ВЮС81аЬ^М1Ы1(0) объясняется как их достаточной эффективностью, так и желанием сопоставить результаты наших тестовых расчетов и аналогичных результатов из работы [4], в которой использовались те же компоненты. Как отмечалось во введении, схема в [4] строилась на основе аппроксимации методом МКЭ уравнений (1), (2) и динамическом исключении одной из неизвестных на шаге метода Ныотона.

0

> -5

100 20 40 60 80 100

Рис. 3. Решение задачи при t = 1 день на разных сетках (nz = 41 и nz = 201)

6. Результаты вычислений

Приводом результаты расчетов четырех тестовых задач (IBM PC. 2.2 GHz. С--) из [4]. сопоставляя результаты вычислений там. где это возможно. Данные из [4] будем указывать в фигурных скобках.

Во всех вычислениях полагалось s* = 1, в = 2, Z = 0.85, 5 = 10-4. Схема МКЭ строилась па основе билинейных конечных элементов на ортогональной сетке из nx х nz узлов (линейных элементов на равномерной сетке из nz узлов в одномерном случае). В каждой зоне однородности среды сетка выбиралась равномерной по каждому направлению с числом узлов, пропорциональным размеру зоны в данном направлении.

1. Однородная одномерная задача. Рассматривается изначально достаточно сухая одномерная колонна длины L = 1 м при следующих параметрах VG:

n = 2, а = 3.35 м-1, ф = 0.368, s* = 0.277, K0 = 0.922 • 10-4 м •с-1, ^|z=Q = -0.75 м, ^|z=L = -10 м, = -10 м.

Мы выбрали, как ив [4], ti = 1 с, время вычисл ений T =1 сут и различное

nz

nz = 41, 201{201}, 1001 точек оказалось равным 251, 468(437}, 624, число итераций Ныотона 275, 485(443}, 626. Соответственно, число неудачных шагов было равно 24, 17(6}, 2. Среднее число итераций BiCGstab на слой по времени составило 4.8, 5.0, 8.3 итераций.

В данной задаче ||u||TO = 0.37, тогда как sa = 0.82. Поэтому переменная u всюду совпадает с se, и метод фактически равносилен аналогичному методу решения s

nz = 1001 nz = 201

дают с графической точностью. Можно отметить хорошую точность решения на грубых сетках. В данном случае метод из [4] оказался на 10% экономичнее по числу итераций Ныотона. Возможно, это объясняется тем, что в [4] использовался метод прогонки для решения системы. Тем не менее, как нам представляется, экономичность реализации шага Ныотона в целом компенсирует эти потери.

Табл. 1

Параметры пористой среды для задачи 2

среда Ко [м/с] [1] а [1/м] п [1]

песок 6.262 • 1СГЬ 0.3658 0.07818 2.80 2.2390

глипа 1.516 • 1СГ6 0.4686 0.2262 1.04 1.3954

Табл. 2

Параметры пористой среды для задачи 3 и 4

зона Ко [м/с] W а [1/м] п [1]

1 9.153 • 10" -ь 0.3680 0.2771 3.34 1.982

2 5.445 • 10" -Б 0.3510 0.2806 3.63 1.632

3 4.805 • 10" -Б 0.3250 0.2643 3.45 1.573

4 4.805 • 10" -4 0.3250 0.2643 3.45 1.573

2. Рассмотрим задачу, в которой реализуется весь диапазон насыщенности. Геометрия задачи определена на рис. 4. параметры неоднородной пористой среды приведены в табл. 1. Граница области непроницаема, кроме участка, на которой задан достаточно большой поток q° = 0.5 м/сут.

T=

= 1 сут. Было выбрано Tmax = 1200 с, т1 = 10{0.9} с. Начальная насыщенность очень мала: = —500 м. В момент времени t = T в средней части области течения образуется область полного насыщения. Приведем результаты вычислений на двух сетках.

х х х

потребовалось 526 {1211} (775) шагов, при этом 41 {345} (54) шагов были отвергнуты. Средний шаг по времени составил 164 (111) с. дисбаланс массы жидкости оказался равным 0.04 {0(10-2}} (0.05)%. Среднее число итераций BiCGstab на слой по времени составило 5.4 (11.0) итерации. Вычисления потребовали 16 с (4 мин 17 с) процессорного времени. Изолинии насыщенности изображены на рис. 6.

В данном расчете наш метод неожиданно оказался существенно экономичнее.

3. Приведем результаты решения другой задачи, область течения которой представлена на рис. 5. а параметры неоднородной пористой среды в табл. 2. Граница области непроницаема, кроме участка, на которой задан слабый поток q° = 0.02 м/сут. Параметр нелинейпости n слабо меняется в области.

Время расчета T = 30 сут, Tmax = 20 мин, т1 = 10{86} с. Начальная насыщенность достаточно мала: = —100 м. Приведем результаты вычислений на двух сетках.

х х х

потребовалось 263 {279} (521) шагов, при этом 1 {0} (28) шаг был отвергнут. Средний шаг по времени составил 0.11 (0.058) сут, дисбаланс массы жидкости оказался равным 0.002 {0(10-3)} (0.005)%. Среднее число итераций BiCGstab на слой по времени составило 3.3 (4.5) итерации. Вычисления потребовали 4 с (87 с) процессорного времени. Изолинии насыщенности изображены на рис. 6.

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

n

n

Яп=0.5м/д

щщщщ

песок

<//о=-500м глина

песок

1 м I

^ 4м

5 м

Рис. 4. Область течения в задаче 2

Рис. 6. Изолинии насыщенности при £ = = 1 сут на разных сетках (51 х 61 - точки, 121 х 121 - сплошные)

I 2 35 ш I

ЛШШШШ1 Ял- 0.02м/д

зона 1

зона 2

зона 4 Д зона 3 У>=-100м

1 к I 2 м

1 1 1 ' 1 1

Ь-^-Ч

Рис. 5. Область течения в задачах 3, 4

Рис. 7. Изолинии насыщенности. Точки -х

х

течения. Как следствие, в момент времени Ь = 30 сут в нижней части области образуется «язык» зоны полного насыщения.

Для интегрирования по времени на сетке из 93x22 {90x21} 120x 120 узлов потребовалось 2191 {1202} 2647 шагов, при этом 0 {13} 30 шагов были отвергнуты. Средний шаг по времени составил 1183-1065 с, дисбаланс массы жидкости оказался равным 0.02 {0(10_3)} 0.03%, среднее число итераций ВЮСз1аЬ на слой по времени составило 2.4 (3.8) итерации. Вычисления потребовали 34 с (6 мин 56 с) процессорного времени. Изолинии насыщенности представлены на рис. 8. В данном тесте наш метод потребовал на 80% больше итераций Ньютона, чем метод из [4]. Мы объясняем это достаточно большим перепадом показателя п в зонах.

Заключение

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

Рис. 8. Линии равной насыщенности при T = 30 сут. Слева — наши вычисления (120 х 120=14400 узлов сетки), справа - из [4] (28917 узлов сетки)

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

Авторы выражают благодарность А.Г. Егорову и М.Ф. Павловой за плодотворные дискуссии, способствовавшие появлению работы.

Работа выполнена при финансовой поддержке РФФИ (проект Х- 07-07-00183).

Summary

A.A. Akhtareev, R.Z. Dautov. Mixed variable technique for simulating uiisat.urat.ed-sat.urat.od flows.

Mixed variable technique for simulating unsat.urat.od-sat.urat.ed flows is proposed. It. based on a new form of governing equation in term of mixed variable and finite elements method with numerical integration. For homogeneous medium new variable coincides with saturation in one regions and with "normalized" pressure head otherwise. Efficient, algorithm for solving discrete problem is proposed. Computational accuracy and robustness of the method are demonstrated by solving four test, problems.

Список литературы

1. Brooks R.H., Corey A.T. Properties of porous media affecting fluid flow // J. Irrig. Drain. Div. Am. Soc. Civil Eng. 1966. V. 92. P. 61 88.

2. Van Genuchten M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils // Soil Sci. Soc. Am. J. 1980. V. 44. P. 892 898.

3. Bixler N.E. An improved time integrator for finite element, analysis // Comm. Appl. Num. Moths. 1989. V. 5. P. 69 78.

4. Diersch, H.-J.G., Perrochet P. On the primary variable switching technique for simulating uiisat.urated-saturat.ed flows // Adv. Water Res. 1999. No 23. P. 271 301.

5. Diersch, H.-J. Finite element, modelling of recirculating density driven saltwater intrusion processes in groundwater // Adv. Water Res. 1988. No 11. P. 25 43.

6. Van der Vorst H.A. Bi-CGSTAB: A fast, and smoothly convergent variant of BiCG for the solution of nonsymmetric linear systems // SIAM J. Sei. St.at. Comp. 1992. No 13. P. 631 644.

7. Meijerink J., Van der Vorst H.A. Guidelines for the usage of incomplete decompositions in solving sets of linear equations as they occur in practical problems // J. Comput. Phys. 1981. No 44. P. 134 155.

Поступила в редакцию 16.10.07

Ахтареев Айдар А. младший научный сотрудник НИИ математики и механики им. Н.Г. Чеботарева Казанского государственного университета.

Даутов Рафаил Замилович доктор физико-математических паук, профессор кафедры вычислительной математики Казанского государственного университета. E-mail: E-mail Rafail.DautovQksu.ru