К.С. Басниев, А.М. Власов, И.Н. Кочина, В.М. Максимов - Подземная гидравлика (1132331), страница 54
Текст из файла (страница 54)
Непосредственная проверка сходимости разностных схем вызывает большие затруднения. В теории разностных схем доказывается, что схема, которая аппроксимирует исходную задачу (когда погрешность аппроксимации стремится к нулю, если стремится к нулю шаг дискретизации) и устойчива (т. е. малым возмущениям начальных данных и разностного оператора соответствуют малые отклонения решений), является сходящейся. Исследования аппроксимации и устойчивости оказываются относительно более простыми. В соответствующих разделах теории разностных схем они описаны достаточно подробно.
Трудности решения системы конечно-разностных уравнений в первую очередь обусловлены ее большой размерностью, равной числу дискретных точек, в которых ищутся значения функций. Размерность системы можно уменьшить увеличением шага дискретизации. Однако этот путь приводит к потере точности, и им нельзя пользоваться без достаточного обоснования. При фиксированной размерности трудоемкость решения системы конечно-разностных уравнений зависит от типа разностной схемы Поясним это на примере системы уравнений (13.3), где во втором 285 слагаемом вместо индекса, указывающего номер временного шага, пока проставлены звездочки.
Рассмотрим два варианта: и+1 и ае — и л и — — (и1+1 — 2ие +и1 1) =О; т аг л+1 л — (и,"е1' — 2иг ' + и,"+1 ) = О, (13,9) т йг 1 = 1, 2, М вЂ” 1; п=1,2 Совокупности узлов, участвующих в разностных схемах, можно изобразить схематично в виде так называемых шаблонов. Для первой и второй схем соответственно шаблоны имеют вид, представленный на рис. 13.5. Ргга, !ЗХ шаблоны для явной (а) и неявной (б) елен ъ-1 в эг! ъ-( в йгв( Первая схема называется явной.
Как наглядно видно из шаблона, в этой схеме каждое значение искомой функции на (п + 1)-м временном слое индивидуально определяется через три значения функции на предыдущем л-м слое. Разрешив (13.8) относительно и,"+', получим л+1 и , 'г ; и и л и, и; + — (иг ~1 — 2ие+ие 1). ьг Используя начальные условия, т. е. значения функции на нулевом слое, мы довольно просто можем вычислить последовательно значения функции на любом последующем временном слое.
Но явная схема имеет очень существенный недостаток: она оказывается сходящейся только при соблюдении ограничительного условия т1йг ( 11'2, т. е. условно устойчива. Вторая схема называется неявной. Она сходится при любом отношении шагов. Но, как видно из шаблона, непосредственно вычислить индивидуальное значение функции на (п + 1)-м слое по известным значениям на и-м слое нельзя. На каждом слое нужно полностью решать всю систему уравнений.
Как уже подчеркивалось ранее, система конечно-разностных уравнений является алгебраической, и поэтому к ней применимы известные методы решения алгебраических уравнений. В то же время отметим, что каждое неявное конечно-разностное уравнение содержит только три значения искомой функции в соседних узлах. Вследствие этого матрица коэффициентов системы конечно-раз- 286 постных уравнений имеет специальный, так называемый трехдиа- гональный вид. Для системы (13.9) матрицей является а Ь 0 Ь и Ь О 0 Ь а Ь 0 0 Ь а Ь 0 Ь а! где а =! + 2т/6', Ь = — т/Ь', а все злементы, не лежащие на трех выделенных диагоналях, равны нулю. Для решения систем уравнений подобного вида разработаны специальные методы, среди которых одним из наиболее популярных является метод прогонки. Рассмотрим его на примере решения системы (13.9) с начальными. н граничными условиями (!3.4), (13.5), (13.7), Для удобства перепишем уравнение в развернутом виде и!+11 — (2+ з) и,"+'+ и1ф' = — и11, з =- Ь'/х.
(13.107 Будем считать, что значения функции в соседних узлах связаны между собой ссютношениямн л+! л+! Г, +1, л!.1! и! = а! (Ь! + и!+1); л+1 л+1 Глл-!.1, л.1.1) и; !=а! ! (Ь! — ! +и1 Подставив (13.12) в (13.10), получим а";+!' (Ь!+!' +и!+') — (2+ а) и!+'+ и,"~+!' = — зи";, откуда следует (13. 13) л+1 Сравнивая (13.13) с (13.11), получаем соотношения «-!-1 ! лл-1-1 л-!-1 «л -!-1, л а! =-, Ь; = — а! '1Ь; ! +ж1, (13. 147 которые по известным значениям а! 1, Ь; '! позволяют вычислить «+1 «+1 и1, ! л+1 Ьл+1 Для узла ! = 1 уравнение (!3.10) и соотношения (13.11) имеют вид и,+ — (2+з)и",+ +из = — ш1, (13.!А (13.16) Используя граничное условие (13.5) из (13.15) получим и1+ — [(лР1~ +и!1) +из+ 1 (13.
17) 2+« 287 и далее в результате сопоставления равенств (13,17) и (13.16) имеем »+ ! , Ь»+! ,»+ + (13.18) 2 -!- » Соотношения ( 1 3. 1 8) и ( 1 3 . 1 4) позволяют последовательно рассчитать пары коэффициентов (а!"+, Ь|+!), (а~~, Ьр~ ), (ам — ! Ьм+'!). Эта процедура известна как «прямой ход». Используя найденные значения пар коэффициентов и второе из граничных условий (13.5), с помощью (13,11) можем вычислить искомые значения функции во всех узлах, начиная с ! =- М вЂ” 1: ам 1 =-ам ! (Ьм — !+»р2 ); им — 2 =ам — 2(Ьм .»+ам !), В этом состоит операция «обратный ход» Повторяя «прямой» и «обратный> ходы на каждом последующем временном шаге, найдем численное решение задачи в любом заданном интервале изменения времени. Метод прогонки удобен тем, что требует относительно небольших объемов оперативной памяти и затрат времени на проведение расчетов.
% З. ЧИСЛЕННОЕ РЕШЕНИЕ ДВУМЕ~НОй ЗАДАЧИ НЕСТАЦИОНАРНОЙ ФИЛЬТРАЦИИ УПРУГОЙ ЖИДКОСТИ В НЕОДНОРОДНОМ ПЛАСТЕ Пусть в горизонтальной плоскости (х, у) имеется область Р(, за- нятая нефтью и содержащая скважины — точечные источники или. стоки. Будем считать, что пласт — неоднородный по проницае- мости: й« =- й» (х, у), а разработка залежи ведется при упругом режиме фильтрации, Для простоты будем предполагать, что область фильтрации Р, имеет форму прямоугольника: Х, < х ( Х„ У, -С у < У, (рис. !3.6). На границах области фильтрации х = Хы х = Х, и у = у» задано соответственно распределение давлений р=р,(у, 1), р=р,(у, 1), р=р,(х, 1). Граница пласта у = У! считается непроницаемой, т, е. на этой границе нормальная составляющая скорости фильтрации (или др/ду) равна нулю.
Пусть в начальный момент времени 1, в пласте (область Р!) задано распределение давления по координатам, т. е. р =- р»(х, у) при 1 =- 1», Тогда задача о нахождении распределения давления р (х, у, () в процессе эксплуатации залежи сводится к решению (интегриро- ванию) дифференциального уравнения параболического типа (типа теплопроводности), с переменными коэффициентами (см. гл. 3, 6), которое можно представить в следующем обобщенном виде д! дк ~ дх/ ду Х ду/ Ь = Ь (х, у), й = й (х, у), 1 =1 (х, у, 1) гвв У1 '=0 2-11'2 ' '+ /2 Рис. И.б.
Схема области фильтрации Рис. И.7. Дискретный аналог не- упругой жидкости прерывной области фильтрации в области Р = Р, Х Рт, Рт — — (1 ) 1,( при следующих начальных и граничных условиях: Рга (хь уг(, хг —.. Й, у; =-1й, 1=-0, тИ, 1=-0, У. Построим, даяее, конечно-разностный аналог уравнения (13.19), используя интегро-интерполяционный метод. Выделим в области Р1 квадрат с центром в точке (х,, йу) и сторонами, образованными отрезками линий х = х, +. 612, д = = — у;~ Ы2 (см. рнс. 13.7).
Рассмотрим тройной интеграл от обеих частей уравнения (13.19): ~а+2 1+п2 /+п2 к у. дн 1у,~х,~( ~ ~ ~ ~д йдо+ к у. 1 — 12 ! — 112 + д и др +~),( ~„,(( др ду га-1-2 1ч-ь2 1+1,2 к. у к у. 1 — 12 1 — 12 289 р = ф (х, у) при (13.20) Р=Ф,(у, 1) при х=Х„ (13.21) р='Фе(у, 1) при к=Хе, (13.22) =0 при у=ум (13.23) р=.ф(х, 1) при у=)ке. (13.24) Здесь искомая функция р (х, у, 1) соответствует давлению, й = йе (х, у)7р; д = ра = тй + (1с — коэффициент упругоемкости пласта; 1 — плотность источников и стоков, моделирующих работу добывающих и нагнетательных скважин.
Со способами моделирования скважин можно познакомиться в работе (12(, Будем реп)ать задачу приближенно с использованием метода конечных разностей. Для этого заменим непрерывную область Р, ее дискретным аналогом — квадратной сеточной областью (рис, 13.7): Выполняя интегрирование по каждому слагаемому в порядке, соответствующем типу производной, получим "(+12"!.=! 2 !л+! !+1(2 Ь~Р + — рл),4уДх = ~ ~ ~(й '1' ) к 2 — 12 ! — (г Ц. л Г (л+! к!+1!2 — ~ — "') 1~ 4(~- ~ ~ ~(й — "') ~К вЂ” — 'Р) 1Ы(-, к ( — ! '2 !и +! «3-(-!(2 9(+1 2 + 5 ) ) Иуухи. 1 — 12 ! — 1(2 Это соотношение — точное. Используя формулы приближенного интегрирования, представим его в следующем виде: [о (Рл ( 1 — />л)) (!> охоу — ~(й — ) — ~/~ — ) 2 2.
(1> ! ! — (( ! х йу~(+ ~(й ~~ ) ' — (й ~~ ) 1 йхЛ(+ (2> л +~(х(, у(, г, ) Лх(куМ, (а>1 Х( ( (Х! ( ( Х(+! (р>! У> ~ (У!' ( ~ У((-(( Г» ~~ ((л ) ~ ~(~.~.2 (а, р, у>=1, 2, 3, 4. Произведения ((др!дх и ((др!ду в точках с полуцелыми индексами заменим дискретными аналогами: (>2 — ) Ж й!+1 г дР ) й Ргь! — Р( (>2 ) л!' — 1 2 где й!2((2 =2й(й~~ l(й!+ >2!л() й( 1 =2й.й,. Дй,.-~-й,г() Подставим полученные выражения в (13.25), предварительно разделив все слагаемые на Лх/гуМ и положив приближенно х!/"' —— =хо у; '=у/, /л'=/и, т.
е. отнеся все средние величины в интегралах к узлу х„у;, /и. В результате получим конечноразностный аналог двумерного уравнения (13,19): и-1-1 л л+1 л+1 л+1 «+1 Рс! Р/,1' Р1ч ь 1' Рс/ Р/,1 Р/,/ 5!, 1 = /21+1 2. 1 «1 — 1!2.! + т аг ' ' 62 и+1 л+! «+! л+1 Рс /+1 — Р1, /, Р/, / Р1, / — ! ил+1 (13.26) 1=1„М вЂ” 1, /=1, Л~ — 1.
Дискретные аналоги начального и граничных условий строятся по ранее рассмотренным схемам и принимают вид при п=О рс/=-1р!,;, 1'=О, М, /=О, /!/; (!3.27) при 1=0 рг,/= — Ф1,/, 1=1, /(/ — 1, п=1, 2, (13.28) при '=М Рм /=% / /=1, /!/ — 1, п=1, 2, (13.29) п=1, 2, (13.30) и = 1 2 и л при !=О рс 1=р,ь !=-1, М вЂ” 1 л л /!И' / Ь / — 1, 2«+1!2 29! при /':= /!/ р/, л = !р/, ! = 1М вЂ” 1, (13.31) Таким путем вместо исходной краевой задачи (13.19) — (13.24) мы получаем конечно-разностную задачу (13.26) — (13.31).
Для решения алгебраической системы уравнений (13.26)— (13.31) можно использовать различные общие и специальные методы. Из числа последних большое распространение получил метод смены направлений, суть которого заключается в следующем. Шаг по времени Л/ =- /„./! — /и разбивается на два половинных шага /«.2! /л+1/2 /л+1/2 /и М!2. На каждом полушаге вместо системы (13.26) — (13.31) решается ее модификация, явная по одному направлению и неявная по другому (направления чередуются). Решаемые системы имеют вид, например, на первом полушаге л+12 л Рл!-112 л+1!2 Р/+1.
/ Р/. / 0,52 /е и+Пг и-~-Пг л и Р/ / — Р; !,'/ „Р/4+1 — Рц — + лп /+1!2 при 1=0 рл! = — Ф! ! — — Л Ф1,ь 1=1, А! — 1, л+1,2 1 + т 2 ' 4 при 1'=--М !гм,! — — — Фг,! — — Л,Ф2,1, 1=1, Л' — 1, л+!!г 1 + т 2 ' 4 где ф-1- 1 (фл-1-1 ( фл). 2 ф — фл+1 Фл — ! . Ф! — !в Лгф! — йь !+1,'2 ~~ — йь ! !'2 Лг ' Лг на втором полушаге л+! «+1!г ,л+!!г .„л+!!г ! ! ! ! й Р!+1,! Р!.! О Зт 1+112 ! Лг л+Ь2 л+1!2 л-! Ы-1 Р1-1,'! Рь !+! — Р1, ! ~! — 1,'2.