— матрица Якоби вектор-функции 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.