Прогнозирование поведения интерфейса при двухфазной фильтрации в реальной области Текст научной статьи по специальности «Математика»

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Справочник по цифровой вычислительной технике / Под ред. Малиновского Б. Н. - Киев: Техника, 1974.-511 с.

2. Чернов Н.И. Линейный синтез цифровых структур АСОиУ. - Таганрог: Изд-во ТРТУ, 2004.-118 с.

3. Искусство схемотехники. Т1. П. Хорвиц, У Хилл. / Под ред. Гальперина М.В.

- М:. Мир, 1983.- 598 с.

В.В.Ершов

ПРОГНОЗИРОВАНИЕ ПОВЕДЕНИЯ ИНТЕРФЕЙСА ПРИ ДВУХФАЗНОЙ ФИЛЬТРАЦИИ В РЕАЛЬНОЙ ОБЛАСТИ

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

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

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

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

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

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

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

др 3 д ( ч

п^~ + Ъ^~(иаР) = 1 ,

д а=1 дха

3 ЯГ 3 я ( дС ^

ді а=і дха а=і дха

Б

V дха у

рдыа др и

— Т и а — ЕР^Х3

п дt дха к

Здесь р(х^) и и(х^) - давление и скорость жидкости, р(С) и /,(С) - ее плотность и вязкость, п(х) и к (х) - пористость и проницаемость среды, С (х^ )

- концентрация одного флюида в другом, Б(х) - коэффициент диффузии, а функции 1(х^) и 8(х^) моделируют внешние источники (осадки, испарение, транспирацию влаги растениями, откачку воды в скважинах, засоление, примеси).

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

Полная модель явления должна учитывать положение свободной поверхности. Считая, что уравнение свободной поверхности имеет вид Е(х^) = х^ — h(Xl,X2,t) = 0, запишем условие на ней в виде

дЕ дЕ дЕ дЕ дh дh дh

------= и 1------Ъ и2-----Ъ из---- или П = из — и 1-----------------и2-----.

дt дх1 дх2 дх3 дt дх1 дх2

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

иа( х^0 ) = иа0(х), Р(х^0) = P0(x), С(х,0) = С0(х1 Мх1,х2^0) = Ых1,х2)-

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

на непроницаемых границах равны нулю нормальные компоненты скорости). Концентрация (или ее нормальная производная) также задана на границе как функция времени.

Рассмотрим область £2^ Я . Введем в £2х [0,Т] равномерные прямоугольные сетки по пространственным переменным и по времени:

®а = Iха I ха = *а^^а , *а = 0, ^а }> = \п \ ^п = п т , п = 0, ^ т }.

Уравнению конвекции-диффузии поставим в соответствие следующую аддитивную разностную схему, аппроксимирующую его в суммарном смысле [2,3]:

ёС(1)

йі

йС(2)

йі

+ К(і)С(1) + БС = 8, і) < і < і)+1, + К(і)С(1) + БС(2) = 8, і) < і < і) +1

с начальными условиями

С<1>(0) = Со, С<2>(іі) = С(1>(1) +1).

Здесь

Б = Е Б

(а)

а=1

о<а)>.=—іа<а>ух

: )

а 'х ха

х є® .

к = Е к

(а)

а=1

К(а)у — ^Ь+ Уха + Ь— аУха , Ь-а < ха < Іа ^а ,

0

ха 0, ха Іа

Дискретизация по времени приводит к чисто неявной схеме для нахождения концентрации флюида:

^п+1/2 _

0,5т

£*п+1 _ £*п+1/2

+ К(іп+1/2)Сп+12 + БСп = 8п,

(1)

+ К(іп+1/2)Сп+12 + БСп+1 = 8п,

п+1

0,5т

где 8п = 8(х,іп

+1/2) • На первом полушаге на верхний слой выносится

конвективный перенос, на втором - диффузионный. Можно показать, что полученная схема безусловно устойчива [2].

Уравнению Дарси соответствует схема, позволяющая найти скорость жидкости на новом временном слое явным образом [3]:

Р

„ п+1 „ п

иа иа

= —Р

—- и ',

а

а

— gР^x3.

(2)

т

Наконец, определим поле давлений. Разностная схема, соответствующая уравнению неразрывности (с учетом уравнений (2)), имеет вид

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

Вычисление скоростей жидкости и положения свободной поверхности по явным схемам [3] труда не представляет, поэтому опишем только методы решения систем (1) и (3). Для каждого узла сетки (кроме граничных) уравнения можно записать в виде

а затем свести их в систему с семидиагональной матрицей Ау = Н, дополнив предварительно аппроксимациями условий на границе области. На практике подобные системы можно решать различными блочно-итерационными методами, такими как В80Я, семейство попеременно-треугольных методов, методы переменных направлений, либо методами расщепления по пространственным переменным (локально-одномерные схемы). Поставленная задача была решена модифицированным итерационным методом переменных направлений, который можно использовать для решения как эллиптических, так и параболических уравнений.

Запишем систему Ау = Н в виде (Х1 + Х2 + Х3 + Т)у = Н , где X1, Х2, Х5, являясь трехдиагональными матрицами, учитывают в А вторые производные по х1, х2 и хз соответственно, а Т - производную по времени

(для системы (3) Т = 0). Тогда итерационная схема, соответствующая одной из модификаций метода Писмэна-Рэчфорда, имеет вид:

Здесь Б - диагональная матрица, содержащая сумму диагональных элементов

Х1, Х2 и Хз. На каждой итерации производится расчет неизвестных

классическим методом прогонки попеременно вдоль трех направлений сетки.

Решение эллиптической системы (3) можно искать по упрощенной схеме:

ауі—1 + Ъу)—1 + суі—1 + йу + еуі+1 + ¿у)+1 + gyl+1 = к,

)Б — Х1 — Х2 — 2Хз

\х1 + г = (г — Х2 — Х3 ^ + Н,

\х2 + г (у+1')Б )у** = (г (у+1Ъ — Х3 Ы — Ху +

В литературе существует несколько подходов к оценке итерационных параметров г. Вполне пригоден подход (предложенный Уайнштейном), при котором гт;п оценивается как минимум выражения (на всей сетке)

( -гг2 -гг2 -гг2 ^

ж ж ж

2N] (1 + Ш] ) 2N2 (l + m2 ) 2N3 (1 + m3 )j

где величины Ш], m2 и т3 зависят от шага сетки и коэффициентов уравнения. Значение rmax, как правило, полагают равным 1 + 2, а промежуточные

итерационные параметры образуют геометрическую прогрессию.

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Bear J., Verruijt A. Modeling Groundwater Flow and Pollution. - D.Reidel Publishing Company, 1987.

2. Самарский А.А., Вабищевич П.Н. Аддитивные схемы для задач математической физики. - М.: Наука, 1999.

3. Афонин А.А., Ершов В.В. О динамике грунтовых вод прибрежных водоносных горизонтов // Известия ТРТУ, 2005.