Смекни!
smekni.com

Н. Г. Чурбанова Москва 2011 г (стр. 4 из 6)

— матрица Якоби вектор-функции F(x). Подставив это

в (4.3), получаем явную формулу метода Ньютона

, (4.4)

Эту формулу, требующую обращения матриц на каждой итерации, можно переписать в неявном виде:

. (4.5)

Применение (4.5) предполагает при каждом k = 0,1,2,... решение линейной алгебраической системы

относительно векторной поправки

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

.

К решению таких линейных систем можно привлекать самые разные методы как прямые, так и итерационные в зависимости от размерности n решаемой задачи и специфики матриц Якоби

(например, можно учитывать их симметричность, разреженность и т.п.).

Сравнивая (4.5) с формальным разложением F(x) в ряд Тейлора

,

видим, что последовательность (

) в методе Ньютона получается в результате подмены при каждом k=0,1,2,... нелинейного уравнения F(x) = 0 или, что то же (при достаточной гладкости F(x)), уравнения

линейным уравнением

т. е. с пошаговой линеаризацией.

3. Алгоритм решения в случае неоднородной среды(двумерная область).

Общий вид расчетной области представлен на рис.2:

Рис.6

В верхней части области находиться источник. А также в области присутствует среда с меньшей проницаемостью.

(3.5)

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

(3.6)

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

3.1 Интерфейсные условия

Рис. 7. Расчетная ячейка на границе раздела сред

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

(3.7)

где

– функция, обратная функции Брукса – Кори в области G2, выражающая зависимость насыщенности от капиллярного давления;
и
– значения водонасыщенности на интерфейсе соответственно со стороны областей G1 и G2;
– водонасыщенность, при которой достигается пороговое давление области G2, то есть она определяется из условия
(см. рис. 8). Таким образом, чтобы DNAPL попал в область G2, его насыщенность в области G1 на интерфейсе должна быть больше, чем
.

(3.8)

=
(3.9)

Рис. 8. Непрерывность капиллярного давления и разрыв насыщенности на интерфейсе.

3.2 Разбиение расчетной области на группы узлов

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

Рис.9 Разбиение расчетной области на группы узлов (расчетные узлы находятся в серединах ячеек)

1-Область расчетные узлы, которой и смежные с ними узлы, полностью лежат в более проницаемой среде -

.

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

.

3-Область расчетные узлы, которой принадлежат верхней границе раздела сред

.

4-Область расчетные узлы, которой принадлежат левой границе раздела сред-

.

5-Область расчетные узлы, которой принадлежат нижней границе раздела сред-

.

6-Область расчетные узлы, которой принадлежат правой границе раздела сред-

.

7-Расчетный узел, находящийся на левом верхнем угле раздела сред-

.

8- Расчетный узел, находящийся на правом верхнем угле раздела сред-

.

9- Расчетный узел, находящийся на левом нижнем угле раздела сред-

.

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

.

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

.

12- Область расчетные узлы, которой и смежные с ними узлы, полностью лежат в менее проницаемой среде-

.

1) Счет в первой и двенадцатой области

Так как все точки разностной схемы для расчета

для точек первой и двенадцатой области принадлежат одной среде, то в них расчет ведется аналогично счету в однородной области (см. §3 пункт 1).

2) Счет во второй и одиннадцатых областях

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

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

Приведем алгоритм расчета для второй области(для одиннадцатой расчет аналогичен, заменяются только параметры среды )

Уравнение неразрывности примет вид:

(3.10)

Тогда разностное уравнение для расчетного узла примет вид:

(3.11)

Скорости фильтрации в полуцелых точках вычисляются по формулам:

(3.12)

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

(3.13)

Далее, выражая

, подставляем получившиеся значения в уравнения (3.4) и приходим, как и в §3 пункте 1, к системе уравнений, которую можно решать методом Ньютона. Далее алгоритм аналогичен описанному в §3 пункте 1.