Конечные разности: от уравнений к вычислениям

Кратко о главном

Разность заменяет производную с погрешностью, которую показывает ряд Тейлора. У разности вперёд и назад первый порядок: вдвое меньший шаг — вдвое меньшая погрешность. У центральной разности и трёхточечной второй разности — второй: погрешность уменьшается вчетверо.
Погрешность аппроксимации — невязка, которую даёт точное решение, подставленное в разностное уравнение. Схема аппроксимирует уравнение (согласована с ним), если невязка стремится к нулю при измельчении сетки и шага.
Устойчивость — свойство не усиливать ошибки от шага к шагу. Аппроксимация вместе с устойчивостью дают сходимость, и погрешность решения убывает с тем же порядком, что и погрешность аппроксимации.
Анализ по фон Нейману раскладывает ошибку на гармоники и находит множитель перехода $G$ каждой за шаг. Схема устойчива, если $|G| \le 1$ для всех гармоник; для явной схемы уравнения давления отсюда $r = \eta\Delta t/\Delta x^2 \le 1/2$ .
Неявная схема (первый порядок по времени) устойчива при любом шаге и не колеблется. Схема Кранка — Николсон тоже устойчива и точнее — второй порядок, — но при крупном шаге быстрые составляющие у неё почти не затухают и меняют знак: после резких изменений решение «звенит».
Исследование сходимости разделяет погрешности сетки и шага. Экстраполяция по сеткам при неизменном шаге сходится к ответу этого шага, а не к точному решению, поэтому шаг проверяют отдельно.

Ключевые соотношения

\frac{p_{i+1} - p_{i-1}}{2\Delta x} = p'_i + \frac{\Delta x^2}{6}\,p'''_i + O(\Delta x^4)

центральная разность: второй порядок

\frac{p_{i+1} - 2p_i + p_{i-1}}{\Delta x^2} = p''_i + \frac{\Delta x^2}{12}\,p^{(4)}_i + O(\Delta x^4)

трёхточечная вторая разность: второй порядок

\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \eta\left[\sigma\,\delta^2 p_i^{n+1} + (1 - \sigma)\,\delta^2 p_i^{n}\right]

схема с весами; $\delta^2 p_i = (p_{i+1} - 2p_i + p_{i-1})/\Delta x^2$ ; $\sigma$ = 0 — явная, 1 — неявная, ½ — Кранка — Николсон

G(\xi) = \frac{1 - 4(1 - \sigma)\,r\sin^2(\xi/2)}{1 + 4\sigma r\sin^2(\xi/2)},\qquad r = \frac{\eta\,\Delta t}{\Delta x^2}

множитель перехода гармоники с фазой $\xi$ на ячейку; устойчивость — $|G| \le 1$

e \approx C_x\,\Delta x^{p} + C_t\,\Delta t^{q}

модель погрешности для исследования сходимости: $p$ и $q$ — порядки по пространству и по времени

Где чаще всего ошибаются

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

После урока вы сможете

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

От производных к разностям

В гидродинамической модели давление известно не везде, а только в отдельных точках — центрах ячеек — и только в отдельные моменты времени: на сетке по пространству и по времени. Обозначим $p_i^n$ давление в точке $x_i$ в момент $t^n$ ; соседние точки отстоят на шаг $\Delta x$ , соседние моменты — на шаг $\Delta t$ . Производных на сетке нет, есть только числа в узлах. Метод конечных разностей заменяет каждую производную комбинацией этих чисел — конечной разностью (finite difference). Набор узлов, которые входят в разностную формулу, называют шаблоном (stencil).

В уроке Уравнение давления для однофазной фильтрации уравнения ячеек получены другим путём — из баланса массы в ячейке. Для одинаковых ячеек однородного пласта оба пути приводят к одному результату: коэффициент связи, делённый на упругую ёмкость ячейки, равен пьезопроводности, делённой на квадрат размера ячейки. Поэтому всё, что теория разностных схем говорит о погрешности и устойчивости, относится и к уравнениям ячеек симулятора.

Уравнение ячейки как разностная схема

S\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = a\left(p_{i+1} - 2p_i + p_{i-1}\right)\quad\Longleftrightarrow\quad \frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \eta\,\frac{p_{i+1} - 2p_i + p_{i-1}}{\Delta x^2}

где:

$S = V_p c_t$ упругая ёмкость ячейки, $V_p = A\,\Delta x\,\phi$ , м³/МПа
$a = 0{,}0853\,kA/(\mu\Delta x)$ коэффициент связи двух одинаковых соседних ячеек, м³/(сут·МПа)
$\eta = 0{,}0853\,k/(\phi\mu c_t)$ пьезопроводность; $\eta/\Delta x^2 = a/S$ , м²/сут
$p_i^{n},\ p_i^{n+1}$ давление в ячейке $i$ в начале и в конце шага, МПа
$\Delta x,\ \Delta t,\ A$ размер ячейки, шаг по времени, площадь грани, м, сут, м²

Единицы: метрические промысловые: $k$ в мД, $\mu$ в мПа·с, $c_t$ в 1/МПа, $p$ в МПа, $t$ в сут; 0,0853 — постоянная закона Дарси в этих единицах (Проницаемость горных пород)

Справедливо, если

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

Справа — центральная вторая разность, умноженная на пьезопроводность: уравнение ячейки и есть разностная схема для уравнения пьезопроводности. Для неоднородного пласта и разных ячеек коэффициенты меняются от ячейки к ячейке, но устройство то же.

Аналогия: скорость по километровым столбам

Пассажир засекает время у каждого километрового столба и делит километр на время между столбами. Получается средняя скорость на километре, но приписать её можно разным моментам. Если отнести её к столбу в начале километра, это «разность вперёд»: при разгоне она завышает скорость в этот момент. Если отнести к середине километра, ошибки от разгона в первой и второй половине почти уравновешиваются — это центральная разность. Поэтому центральные формулы точнее: их ошибки симметричны и сокращаются.

Ряд Тейлора: погрешность аппроксимации и порядок

Насколько разность отличается от производной, показывает ряд Тейлора (Taylor series): значение гладкой функции в соседней точке выражается через значение и производные в данной точке.

Ряд Тейлора для соседних узлов

p_{i\pm1} = p_i \pm \Delta x\,p'_i + \frac{\Delta x^2}{2}\,p''_i \pm \frac{\Delta x^3}{6}\,p'''_i + \frac{\Delta x^4}{24}\,p^{(4)}_i \pm \ldots

где:

$p_i,\ p_{i\pm1}$ значения функции в узле $i$ и в соседних узлах, МПа
$p'_i,\ p''_i,\ p'''_i,\ p^{(4)}_i$ первая — четвёртая производные по $x$ в узле $i$ , МПа/м, МПа/м², …
$\Delta x$ шаг сетки, м

Единицы: любые согласованные

Справедливо, если

функция гладкая: у неё есть непрерывные производные нужного порядка между узлами;
шаг сетки постоянный (неравномерная сетка — в дополнительном разделе).

Верхние знаки — для узла $i + 1$ , нижние — для $i - 1$ . Чем меньше шаг, тем быстрее убывают члены ряда, и тем меньше отбрасываемая часть.

Вычтем $p_i$ и поделим на шаг — получим разность вперёд; вычтем разложения друг из друга — центральную разность; сложим — вторую разность. Члены, которые остаются после производной, и есть погрешность формулы.

Разностные формулы и их погрешность

\begin{aligned} &\frac{p_{i+1} - p_i}{\Delta x} = p'_i + \frac{\Delta x}{2}\,p''_i + O(\Delta x^2) &&\text{(вперёд)}\\ &\frac{p_i - p_{i-1}}{\Delta x} = p'_i - \frac{\Delta x}{2}\,p''_i + O(\Delta x^2) &&\text{(назад)}\\ &\frac{p_{i+1} - p_{i-1}}{2\Delta x} = p'_i + \frac{\Delta x^2}{6}\,p'''_i + O(\Delta x^4) &&\text{(центральная)}\\ &\frac{p_{i+1} - 2p_i + p_{i-1}}{\Delta x^2} = p''_i + \frac{\Delta x^2}{12}\,p^{(4)}_i + O(\Delta x^4) &&\text{(вторая разность)} \end{aligned}

где:

$O(\Delta x^m)$ члены порядка $\Delta x^m$ : при уменьшении шага вдвое они уменьшаются примерно в $2^m$ раз, —
$p'_i,\ p''_i,\ p'''_i,\ p^{(4)}_i$ производные в узле $i$ , МПа/м, МПа/м², …

Единицы: любые согласованные

Справедливо, если

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

При вычитании разложений сокращаются члены с чётными степенями шага, при сложении — с нечётными. Поэтому у центральной разности и у второй разности погрешность начинается со второй степени шага, а у разностей вперёд и назад — с первой: они ошибаются на величину, пропорциональную шагу и кривизне функции.

Погрешность аппроксимации: (truncation error) — разность между разностной формулой и производной, которую она заменяет, при подстановке точной функции. Здесь будем обозначать её $\psi$ .
Порядок аппроксимации: (order of accuracy) — показатель степени шага в главном члене погрешности: формула имеет порядок $m$ , если $\psi = O(\Delta x^m)$ . При уменьшении шага вдвое погрешность убывает примерно в $2^m$ раз.
Главный член: первый неисчезающий член ряда в погрешности, например $\tfrac{\Delta x}{2}p''$ у разности вперёд. По нему видно не только порядок, но и знак ошибки: при $p'' < 0$ разность вперёд занижает производную.

Разностные формулы для производных

Формула: вперёд, $p'$
Узлы: $i,\ i + 1$
Главный член погрешности: $\tfrac{\Delta x}{2}\,p''$
Порядок: 1

Формула: назад, $p'$
Узлы: $i - 1,\ i$
Главный член погрешности: $-\tfrac{\Delta x}{2}\,p''$
Порядок: 1

Формула: центральная, $p'$
Узлы: $i \pm 1$
Главный член погрешности: $\tfrac{\Delta x^2}{6}\,p'''$
Порядок: 2

Формула: вторая разность, $p''$
Узлы: $i,\ i \pm 1$
Главный член погрешности: $\tfrac{\Delta x^2}{12}\,p^{(4)}$
Порядок: 2

Формула	Узлы	Главный член погрешности	Порядок
вперёд, $p'$	$i,\ i + 1$	$\tfrac{\Delta x}{2}\,p''$	1
назад, $p'$	$i - 1,\ i$	$-\tfrac{\Delta x}{2}\,p''$	1
центральная, $p'$	$i \pm 1$	$\tfrac{\Delta x^2}{6}\,p'''$	2
вторая разность, $p''$	$i,\ i \pm 1$	$\tfrac{\Delta x^2}{12}\,p^{(4)}$	2

Равномерная сетка, гладкая функция. Формулы по времени устроены так же: разность $(p^{n+1} - p^n)/\Delta t$ — это разность вперёд, если её относят к моменту $t^n$ , и центральная — если к середине шага.

Расчётный примерПроизводные давления у скважины по трём точкам

Исходные данные (условные)

Профиль давления (скважина урока Режимы фильтрации): $p(r) - p_{wf} = m\ln(r/r_w)$ ; $m$ = 0,388 МПа = 388 кПа; $r_w$ = 0,1 м
Точка и шаги: $r$ = 10 м; $\Delta r$ = 2 м и 1 м
Точные производные: $p' = m/r$ = 38,8 кПа/м; $p'' = -m/r^2$ = −3,88 кПа/м²

Решение

Разности давления при шаге 2 м: $p(12) - p(10) = 388\ln 1{,}2 = 70{,}741$ кПа; $p(10) - p(8) = 388\ln 1{,}25 = 86{,}580$ кПа.
Разность вперёд: $70{,}741/2 = 35{,}37$ кПа/м, ошибка $-3{,}43$ кПа/м (−8,8 %). Назад: $86{,}580/2 = 43{,}29$ кПа/м, ошибка $+4{,}49$ кПа/м. Центральная: $(70{,}741 + 86{,}580)/4 = 39{,}33$ кПа/м, ошибка $+0{,}530$ кПа/м (+1,37 %).
Шаг 1 м: $p(11) - p(10) = 36{,}980$ кПа, $p(10) - p(9) = 40{,}880$ кПа. Вперёд — 36,98 кПа/м (ошибка $-1{,}82$ ), назад — 40,88 (ошибка $+2{,}08$ ), центральная — 38,93 кПа/м (ошибка $+0{,}130$ , 0,34 %).
Во сколько раз уменьшились ошибки при вдвое меньшем шаге: у разности вперёд $3{,}43/1{,}82 = 1{,}9$ — первый порядок; у центральной $0{,}530/0{,}130 = 4{,}1$ — второй.
Главные члены ряда Тейлора при шаге 2 м: у разности вперёд $\tfrac{\Delta r}{2}p'' = -3{,}88$ кПа/м, у центральной $\tfrac{\Delta r^2}{6}p''' = +0{,}517$ кПа/м (здесь $p''' = 2m/r^3$ ). Они близки к фактическим ошибкам; остаток дают следующие члены ряда.
Вторая производная по трём точкам: при шаге 2 м $(70{,}741 - 86{,}580)/2^2 = -3{,}960$ кПа/м² (ошибка 2,1 %), при шаге 1 м $(36{,}980 - 40{,}880)/1^2 = -3{,}900$ кПа/м² (0,5 %): ошибка тоже уменьшилась примерно вчетверо.

Ответ: При шаге 1 м разность вперёд ошибается на 4,7 %, центральная — на 0,34 %. При уменьшении шага вдвое ошибка первой уменьшилась в 1,9 раза, второй — в 4,1 раза.

Что это значит. Центральная разность с шагом 2 м точнее разности вперёд с шагом 1 м — выигрыш от порядка больше, чем от измельчения. Знак ошибок предсказывает главный член: кривизна профиля отрицательна, поэтому разность вперёд занижает наклон, а назад — завышает. Вблизи скважины производные давления велики и быстро растут к стволу (третья производная — как $1/r^3$ ), поэтому равномерная сетка здесь неэффективна: ячейки у скважины сгущают (Сетки в гидродинамических моделях), а сам перепад у ствола переносят в модель скважины (Модель скважины в ячейке сетки). Данные примера условные.

Что меняется, если…

…шаг уменьшить вдвое — погрешность разности вперёд уменьшится примерно вдвое, центральной — вчетверо, но число узлов вырастет вдвое;
…функция меняется резко (велики высшие производные) — ошибаются все формулы, и помогает только мелкий шаг именно там, где она меняется: у скважин, на фронтах, на границах слоёв;
…значения в узлах содержат погрешность $\varepsilon$ (округление, шум замеров) — к ошибке формулы добавляется ошибка порядка $\varepsilon/\Delta x$ , и при очень малом шаге производная портится. Поэтому производную давления по замерам гидродинамических исследований скважин считают по точкам, разнесённым на интервал сглаживания (Диагностика по производной давления).

Схемы для уравнения давления и их погрешность

Возьмём одномерное линейное уравнение давления для однородного пласта — уравнение пьезопроводности $\partial p/\partial t = \eta\,\partial^2 p/\partial x^2$ (Уравнение пьезопроводности). Вторую производную по $x$ заменяем второй разностью; остаётся решить, давления какого момента в неё подставлять. В уроке Уравнение давления для однофазной фильтрации разобраны два крайних варианта — начало шага (явная схема) и конец шага (неявная). Их удобно объединить в одну схему с весами (weighted scheme): правая часть — взвешенное среднее вторых разностей начала и конца шага.

Схема с весами для уравнения давления

\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \eta\left[\sigma\,\frac{p_{i+1}^{n+1} - 2p_i^{n+1} + p_{i-1}^{n+1}}{\Delta x^2} + (1 - \sigma)\,\frac{p_{i+1}^{n} - 2p_i^{n} + p_{i-1}^{n}}{\Delta x^2}\right]

где:

$\sigma$ вес конца шага: 0 — явная схема, 1 — неявная, ½ — схема Кранка — Николсон, —
$\eta$ пьезопроводность, $0{,}0853\,k/(\phi\mu c_t)$ , м²/сут
$p_i^{n},\ p_i^{n+1}$ давление в узле (ячейке) $i$ в начале и в конце шага, МПа
$\Delta x,\ \Delta t$ шаг сетки и шаг по времени, м, сут

Единицы: метрические промысловые (любые согласованные)

Справедливо, если

линейное уравнение, однородный пласт, равномерная сетка; для неоднородного пласта то же записывают для уравнений ячеек с коэффициентами $a_{ij}$ и $S_i$ ;
при $\sigma > 0$ в каждом уравнении три неизвестных конца шага, и на шаге решают трёхдиагональную систему, как в неявной схеме.

Вес $\sigma$ задаёт, по давлениям какого момента считаются перетоки на шаге. При $\sigma = \tfrac12$ это перетоки по среднему давлению шага — та промежуточная схема, о которой упоминает урок об уравнении давления.

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

Погрешность аппроксимации схемы с весами

\psi_i^{n} = \left(\tfrac12 - \sigma\right)\Delta t\,\frac{\partial^2 p}{\partial t^2} - \frac{\eta\,\Delta x^2}{12}\,\frac{\partial^4 p}{\partial x^4} + O\!\left(\Delta t^2 + \Delta x^4\right)

где:

$\psi_i^{n}$ погрешность аппроксимации: левая часть схемы минус правая при подстановке точного решения, МПа/сут
$\partial^2 p/\partial t^2,\ \partial^4 p/\partial x^4$ производные точного решения в середине шага, МПа/сут², МПа/м⁴
$\sigma,\ \eta,\ \Delta t,\ \Delta x$ вес, пьезопроводность, шаги, —, м²/сут, сут, м

Единицы: любые согласованные

Справедливо, если

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

Член с $\Delta t$ исчезает только при $\sigma = \tfrac12$ . Поэтому явная и неявная схемы имеют первый порядок по времени и второй по пространству, а схема Кранка — Николсон — второй по обоим. Погрешность по пространству у всех трёх одна и та же: её даёт вторая разность.

Три схемы как частные случаи схемы с весами

Схема: явная
$\sigma$: 0
Порядок по $t$: 1
Порядок по $x$: 2
Устойчивость: при $r \le 1/2$

Схема: неявная
$\sigma$: 1
Порядок по $t$: 1
Порядок по $x$: 2
Устойчивость: при любом $r$

Схема: Кранка — Николсон
$\sigma$: ½
Порядок по $t$: 2
Порядок по $x$: 2
Устойчивость: при любом $r$

Схема	$\sigma$	Порядок по $t$	Порядок по $x$	Устойчивость
явная	0	1	2	при $r \le 1/2$
неявная	1	1	2	при любом $r$
Кранка — Николсон	½	2	2	при любом $r$

$r = \eta\Delta t/\Delta x^2$ — сеточное число Фурье. Условия устойчивости выведены ниже анализом по фон Нейману; для неоднородной сетки явной схемы удобнее условие через характерные времена ячеек из урока Уравнение давления для однофазной фильтрации.

Аппроксимация, устойчивость, сходимость

Аппроксимация (согласованность, consistency): погрешность аппроксимации $\psi$ стремится к нулю, когда шаги $\Delta x$ и $\Delta t$ стремятся к нулю. Схема аппроксимирует именно то уравнение, которое нужно решить.
Устойчивость (stability): ошибки, внесённые на любом шаге — в начальных данных, при округлении, при неточном решении уравнений шага, — при дальнейшем счёте не нарастают неограниченно: за конечное время расчёта они остаются меньше своей начальной величины, умноженной на постоянное число, как бы малы ни были шаги (устойчивость разностной схемы).
Сходимость (convergence): при измельчении сетки и шага решение схемы стремится к точному решению уравнения в каждой точке и в каждый момент.

Связь между ними видна из уравнения для ошибки решения $e = p_{\text{точн}} - p$ . Точное решение удовлетворяет схеме с невязкой $\psi$ , расчётное — без неё. Вычтя одно из другого, получим ту же схему для ошибки, но с «источником» $\psi$ на каждом шаге:

Накопление ошибки в устойчивой схеме

\begin{aligned} &e^{n+1} = \mathbf{A}\,e^{n} + \Delta t\,\mathbf{B}\,\psi^{n},\qquad \mathbf{B} = \left(\mathbf{I} - \sigma\Delta t\,\mathbf{L}\right)^{-1},\qquad \mathbf{A} = \mathbf{B}\left(\mathbf{I} + (1 - \sigma)\Delta t\,\mathbf{L}\right)\\ &\left\|e^{n}\right\| \le \left\|e^{0}\right\| + t^{n}\max_k\left\|\psi^{k}\right\| \end{aligned}

где:

$e^{n},\ \psi^{n}$ векторы ошибки решения и погрешности аппроксимации по всем узлам на шаге $n$ , МПа, МПа/сут
$\mathbf{A}$ оператор шага: во что схема превращает ошибку начала шага, —
$\mathbf{L}$ матрица вторых разностей, умноженная на пьезопроводность, 1/сут
$\mathbf{I}$ единичная матрица, —
$\|\cdot\|,\ t^{n}$ норма вектора (например, среднеквадратичная) и время от начала расчёта, —, сут

Единицы: любые согласованные

Справедливо, если

линейная задача; оценка справа верна, если $\|\mathbf{A}\| \le 1$ (схема устойчива) и $\|\mathbf{B}\| \le 1$ — так у схемы с весами для уравнения давления;
при неустойчивой схеме $\|\mathbf{A}\| > 1$ , и вклад каждого шага умножается на $\|\mathbf{A}\|$ на всех последующих шагах.

В устойчивой схеме ошибка за время $t^n$ не больше начальной плюс время, умноженное на наибольшую невязку. Невязка убывает при измельчении (аппроксимация), значит, убывает и ошибка (сходимость) — причём с тем же порядком.

Оценка иллюстрирует главную теорему теории разностных схем: для линейной задачи с корректной постановкой и схемы, которая её аппроксимирует, устойчивость необходима и достаточна для сходимости — теорема эквивалентности Лакса (Lax equivalence theorem). Обе половины важны. Неустойчивая схема не сходится, даже если аппроксимирует уравнение: ошибки округления порядка $10^{-15}$ растут на каждом шаге в фиксированное число раз и рано или поздно заполняют весь результат. Устойчивая схема без аппроксимации может сойтись — но к решению другого уравнения (пример — в дополнительном разделе).

Две разные «сходимости»

Симуляторы сообщают о сходимости итераций на шаге: уравнения ячеек решены с заданной невязкой — сохранение массы выполнено, итерации Ньютона закончились (Схемы IMPES и полностью неявная). Это не сходимость схемы. Идеально решённые разностные уравнения на грубой сетке с крупным шагом всё равно отличаются от решения дифференциального уравнения на погрешность схемы. Её проверяют только измельчением сетки и шага.

Анализ устойчивости по фон Нейману

Как проверить, что $\|\mathbf{A}\| \le 1$ , не вычисляя матриц? Для линейного уравнения с постоянными коэффициентами на равномерной сетке ошибку можно разложить на гармоники — синусоиды разной длины волны — и проследить каждую отдельно. Схема переводит гармонику в ту же гармонику, умноженную на число $G$ — множитель перехода (amplification factor). Метод предложил Джон фон Нейман в 1940-х годах; его называют также методом гармоник или спектральным признаком устойчивости (анализ по фон Нейману).

Пусть ошибка в узле с номером $j$ равна $G^n\cos(j\xi)$ , где $\xi$ — сдвиг фазы гармоники от узла к узлу. При малом $\xi$ гармоника плавная (длинная волна), при $\xi = \pi$ — «шахматная»: соседние значения равны по величине и противоположны по знаку; короче волну на сетке не изобразить. Вторая разность такой гармоники — та же гармоника, умноженная на $-4\sin^2(\xi/2)/\Delta x^2$ (это следует из формулы суммы косинусов). Подставив гармонику в схему с весами, получим:

Множитель перехода схемы с весами

G(\xi) = \frac{1 - 4(1 - \sigma)\,r\sin^2(\xi/2)}{1 + 4\sigma\,r\sin^2(\xi/2)},\qquad r = \frac{\eta\,\Delta t}{\Delta x^2}

где:

$G(\xi)$ во сколько раз амплитуда гармоники меняется за шаг; $G < 0$ — гармоника меняет знак на каждом шаге, —
$\xi$ сдвиг фазы гармоники на одну ячейку, $0 < \xi \le \pi$ ; длина волны $2\pi\Delta x/\xi$ , рад
$r$ сеточное число Фурье: шаг по времени в единицах времени выравнивания давления на одной ячейке, —
$\sigma$ вес схемы, —

Единицы: безразмерная

Справедливо, если

линейное уравнение с постоянными коэффициентами, равномерная сетка; границы не учитываются (бесконечная или периодическая сетка);
для неоднородного пласта условие проверяют по месту — с местными $\eta$ и $\Delta x$ ; для нелинейного уравнения — с коэффициентами, «замороженными» на шаге; тогда условие необходимое, но не всегда достаточное.

Точное решение уменьшает ту же гармонику за шаг в $\e^{-r\xi^2}$ раз. Схема устойчива, если $|G| \le 1$ для всех $\xi$ , и точна для тех гармоник, у которых $G$ близок к $\e^{-r\xi^2}$ .

Условия устойчивости схемы с весами

\sigma = 0:\ \ r \le \tfrac12;\qquad 0 < \sigma < \tfrac12:\ \ r \le \frac{1}{2(1 - 2\sigma)};\qquad \sigma \ge \tfrac12:\ \ \text{при любом } r

где:

$r$ сеточное число Фурье $\eta\Delta t/\Delta x^2$ , —
$\sigma$ вес схемы, —

Единицы: безразмерная

Справедливо, если

получены из $G \ge -1$ для самой короткой, шахматной гармоники ( $\xi = \pi$ ): $G \le 1$ выполняется всегда;
в двух и трёх измерениях для явной схемы $r_x + r_y + r_z \le \tfrac12$ — то же условие, что для равномерной сетки в уроке об уравнении давления.

Для явной схемы $r \le \tfrac12$ означает $\Delta t \le \Delta x^2/(2\eta)$ — предел из урока Уравнение давления для однофазной фильтрации. Схемы с $\sigma \ge \tfrac12$ устойчивы при любом шаге, но это ещё не значит, что они хорошо гасят быстрые гармоники (ниже — о схеме Кранка — Николсон).

Анализ гармоник связывает два взгляда на устойчивость из предыдущих уроков. Каждая гармоника ведёт себя как отдельный «бак» со своим характерным временем $\tau_\xi = \Delta x^2/(4\eta\sin^2(\xi/2))$ : явный шаг умножает её на $1 - \Delta t/\tau_\xi$ , ровно как депрессию одного бака в уроке о гидродинамической модели, и устойчив при $\Delta t \le 2\tau_\xi$ . Самый «быстрый» бак — шахматная гармоника с $\tau_\pi = \Delta x^2/(4\eta)$ ; отсюда $\Delta t \le \Delta x^2/(2\eta)$ — тот же предел, что дал шахматный узор в уроке об уравнении давления.

РасчётМножитель перехода за шаг для гармоник разной длины волны при

r \approx 4

(ячейки 250 м, шаг 30 сут, пьезопроводность 8412 м²/сут) и точное затухание

\e^{-r\xi^2}

. Явная схема выходит за границу устойчивости уже у довольно длинных волн; неявная гасит короткие волны, а длинные — слабее точного; у схемы Кранка — Николсон короткие волны затухают медленно и меняют знак.

Расчётный примерМножители перехода трёх схем на сетке с ячейками 250 м

Исходные данные (условные)

Сетка и порода (урок Режимы фильтрации): $\Delta x$ = 250 м; $\eta$ = 8412 м²/сут
Шаги по времени: 3, 6 и 30 сут
Гармоники: плавная $\xi = \pi/8$ (волна 16 ячеек = 4 км) и шахматная $\xi = \pi$

Решение

Числа $r = \eta\Delta t/\Delta x^2$ : 0,4038; 0,8076; 4,04. Предел явной схемы $r \le 0{,}5$ , то есть $\Delta t \le \Delta x^2/(2\eta) = 3{,}71$ сут.
Плавная гармоника, $\sin^2(\pi/16) = 0{,}0381$ . Шаг 3 сут: явная 0,9385, неявная 0,9421, Кранка — Николсон 0,9404; точно $\e^{-r\xi^2} = 0{,}9396$ — все три близки к точному.
Шахматная гармоника, $\sin^2(\pi/2) = 1$ . Шаг 3 сут: явная $1 - 4\cdot0{,}4038 = -0{,}615$ — затухает, меняя знак; неявная 0,382; Кранка — Николсон 0,106.
Шаг 6 сут: плавная гармоника у явной схемы 0,877 при точных 0,883 — по-прежнему хорошо; шахматная $1 - 4\cdot0{,}8076 = -2{,}23$ : за 10 шагов она вырастает примерно в 3000 раз, за 30 шагов — примерно в 2,8·10¹⁰ раз.
Шаг 30 сут, неявная схема: плавная 0,619 при точных 0,536, шахматная 0,058 — гасится почти полностью.
Шаг 30 сут, схема Кранка — Николсон: плавная 0,530 — почти точно; шахматная $(1 - 2\cdot4{,}04)/(1 + 2\cdot4{,}04) = -0{,}78$ : через 5 шагов от неё остаётся 0,29 начальной амплитуды (знак меняется на каждом шаге), тогда как в точном решении она исчезает за первый же шаг.

Ответ: Явная схема на этой сетке устойчива при шаге до 3,71 сут; при 6 сут шахматная гармоника растёт в 2,23 раза за шаг. При шаге 30 сут неявная схема гасит шахматную гармонику (0,058), а схема Кранка — Николсон почти нет (−0,78), зато плавную передаёт точнее (0,530 против 0,619 при точных 0,536).

Что это значит. Таблица множителей объясняет поведение схем без расчёта по ячейкам. Неустойчивость — это $|G| > 1$ у коротких волн. Запаздывание неявной схемы — это её $G$ для длинных волн, больший точного: возмущение затухает медленнее, чем должно. «Звон» схемы Кранка — Николсон — это $G$ , близкий к −1, у коротких волн. Начальные данные и ошибки округления содержат все гармоники, поэтому неустойчивая гармоника рано или поздно проявится, даже если в первые шаги её не видно. Данные примера условные.

Границы и неоднородность анализ гармоник не учитывает. Обычно они делают предел немного мягче: в полосе урока Уравнение давления для однофазной фильтрации наименьшее характерное время ячеек — 3,3 сут, а точный предел — около 4,1 сут. Поэтому на практике пользуются местным условием по самым мелким и самым проводящим ячейкам: оно надёжно и проверяется для каждой ячейки отдельно.

Схема Кранка — Николсон

При $\sigma = \tfrac12$ перетоки на шаге берут как среднее перетоков по давлениям начала и конца шага — по сути, интеграл по времени считают по правилу трапеций. Это схема Кранка — Николсон (Crank–Nicolson scheme), предложенная в 1947 году для уравнения теплопроводности. В записи уравнения ячейки:

Схема Кранка — Николсон для уравнения ячейки

S_i\,\frac{p_i^{n+1} - p_i^{n}}{\Delta t} = \frac12\sum_j a_{ij}\left(p_j^{n+1} - p_i^{n+1}\right) + \frac12\sum_j a_{ij}\left(p_j^{n} - p_i^{n}\right) - q_i

где:

$S_i,\ a_{ij}$ упругая ёмкость ячейки и коэффициенты связи с соседями, м³/МПа, м³/(сут·МПа)
$q_i$ отбор скважиной из ячейки в пластовых условиях, м³/сут
$p^{n},\ p^{n+1}$ давления начала и конца шага, МПа

Единицы: метрические промысловые, как в уроке об уравнении давления

Справедливо, если

линейное уравнение; отбор $q_i$ постоянен на шаге;
система для давлений конца шага трёхдиагональная (в одном измерении) с диагональю $S_i/\Delta t + \tfrac12\sum_j a_{ij}$ .

Работа на шаге та же, что у неявной схемы, — одна система уравнений, — а погрешность по времени второго порядка. Схема устойчива при любом шаге: $|G| \le 1$ для всех гармоник.

Слабое место видно из множителя перехода: при $\sigma = \tfrac12$ и шахматной гармонике $G = (1 - 2r)/(1 + 2r)$ . При крупном шаге ( $r \gg 1$ ) он стремится к −1: короткие волны почти не затухают и меняют знак на каждом шаге. Пока решение гладкое, коротких волн в нём мало, и это незаметно. Но после резкого изменения — пуска или остановки скважины, смены дебита, резкого начального перепада — короткие волны есть, и у ячеек рядом с местом изменения появляются затухающие пилообразные колебания. Неявная схема такие волны гасит сразу ( $G \to 0$ ): это свойство сильного затухания быстрых составляющих называют L-устойчивостью (L-stability).

Расчётный примерОдин шаг 10 сут по схеме Кранка — Николсон на полосе из четырёх ячеек

Исходные данные (условные)

Полоса (урок Уравнение давления для однофазной фильтрации): $S$ = 316,875 м³/МПа; $a_{12}$ = 68,24, $a_{23}$ = 27,296, $a_{34}$ = 28,433 м³/(сут·МПа); отбор из ячейки 4 — 130 м³/сут в пластовых условиях; начальное давление 20 МПа
Шаг: $\Delta t$ = 10 сут — втрое больше наименьшего характерного времени ячеек (3,3 сут)

Решение

В начале шага давление однородно, перетоков нет, поэтому для изменений давления $\delta_i$ система та же, что у неявной схемы, но с половинными коэффициентами связи: $\tfrac12 a$ = 34,120; 13,648; 14,217. Диагональ $b_i = S/\Delta t + \tfrac12\sum_j a_{ij}$ при $S/\Delta t = 31{,}6875$ : 65,8075; 79,4555; 59,5525; 45,9045. Правая часть $(0;\ 0;\ 0;\ -130)$ .
Прямой ход прогонки: $e_1 = 34{,}120/65{,}8075 = 0{,}5185$ ; $w_2 = 79{,}4555 - 34{,}120\cdot0{,}5185 = 61{,}76$ , $e_2 = 13{,}648/61{,}76 = 0{,}2210$ .
$w_3 = 59{,}5525 - 13{,}648\cdot0{,}2210 = 56{,}54$ , $e_3 = 14{,}217/56{,}54 = 0{,}2515$ ; $w_4 = 45{,}9045 - 14{,}217\cdot0{,}2515 = 42{,}33$ .
Обратный ход: $\delta_4 = -130/42{,}33 = -3{,}071$ ; $\delta_3 = 0{,}2515\cdot\delta_4 = -0{,}772$ ; $\delta_2 = 0{,}2210\cdot\delta_3 = -0{,}171$ ; $\delta_1 = 0{,}5185\cdot\delta_2 = -0{,}089$ МПа.
Давления через 10 сут: 19,91; 19,83; 19,23; 16,93 МПа. Точное решение уравнений ячеек: 19,92; 19,82; 19,11; 17,04 МПа. Неявная схема с тем же шагом давала в ячейке 4 17,39 МПа.
Шаг 20 сут: давление в ячейке 3 по схеме Кранка — Николсон отклоняется от точного попеременно в разные стороны — на +0,24, −0,15 и +0,03 МПа через 20, 40 и 60 сут; у неявной схемы отклонения одного знака — +0,06, +0,11 и +0,08 МПа.

Ответ: Через 10 сут в ячейке со скважиной схема Кранка — Николсон даёт 16,93 МПа при точных 17,04 МПа: ошибка −0,11 МПа против +0,35 МПа у неявной схемы с тем же шагом.

Что это значит. Второй порядок по времени заметен в ячейке со скважиной: при той же работе — одна трёхдиагональная система на шаг — ошибка в ней втрое меньше. В соседней ячейке 3 картина обратная: уже через 10 сут схема Кранка — Николсон ошибается сильнее неявной (+0,12 против −0,06 МПа), а при шаге 20 сут давление в ней «качается» вокруг точного решения, и первые отклонения больше, чем у неявной схемы. Пуск скважины — резкое изменение, а быстрые гармоники у этой схемы затухают медленно и меняют знак. Неявная схема в ячейке со скважиной ошибается сильнее, но её ошибки одного знака, без колебаний, — для прогнозов с частыми изменениями режимов скважин это важнее. Данные примера условные.

РасчётПроверочная задача: в закрытой полосе длиной 2 км давление в левой половине на 2 МПа выше, чем в правой. Профиль через 10 сут, рассчитанный одним шагом 10 сут на сетке из 20 ячеек по 100 м (

r \approx 8{,}4

). Неявная схема даёт гладкий профиль, схема Кранка — Николсон «переворачивает» перепад в двух средних ячейках; на следующем шаге он перевернётся обратно.

Колебания убирают несколькими способами. Надёжнее всего после каждого резкого изменения сделать два-три шага неявной схемой и только потом перейти к $\sigma = \tfrac12$ : неявные шаги гасят короткие волны. Помогает и уменьшение шага после изменений режима. Вес чуть больше половины (например, 0,55) ослабляет «звон» — при крупном шаге шахматная гармоника умножается уже не на −1, а примерно на $-(1 - \sigma)/\sigma \approx -0{,}82$ — и сохраняет малую погрешность по времени, но колебаний полностью не убирает. В гидродинамических симуляторах по времени обычно используют неявную схему первого порядка: многофазные уравнения нелинейны, режимы скважин меняются часто, и надёжность важнее второго порядка, а нужную точность обеспечивают выбором шага (Схемы IMPES и полностью неявная).

Исследование сходимости по сетке и по шагу

Проверку сходимости по сетке — расчёты на нескольких сетках с двукратным измельчением, наблюдаемый порядок и экстраполяцию к бесконечно мелкой сетке — разбирает урок о сетках. Теория разностных схем добавляет к ней две вещи: чего ждать от наблюдаемого порядка и как не перепутать погрешность сетки с погрешностью шага по времени. Для гладкого решения погрешность результата складывается из двух частей:

Модель погрешности для исследования сходимости

e = R - R_{\text{точн}} \approx C_x\,\Delta x^{p} + C_t\,\Delta t^{q}

где:

$R,\ R_{\text{точн}}$ результат расчёта (давление, дебит, время прорыва) и его точное значение, в единицах результата
$p,\ q$ порядки по пространству и по времени: для уравнения давления $p$ = 2; $q$ = 1 у неявной схемы и 2 у схемы Кранка — Николсон, —
$C_x,\ C_t$ коэффициенты, зависящие от решения (его производных) и от схемы, но не от шагов, —

Единицы: в единицах результата

Справедливо, если

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

Если измельчать только сетку, погрешность стремится не к нулю, а к $C_t\Delta t^q$ ; если только шаг — к $C_x\Delta x^p$ . Поэтому шаг и сетку проверяют по отдельности или измельчают согласованно.

Сетка при неизменном шаге. Разности результатов на трёх сетках по-прежнему показывают порядок $p$ — постоянная часть $C_t\Delta t^q$ в них сокращается, — но экстраполяция даёт ответ для этого шага, а не точный.
Шаг на неизменной сетке. Погрешность убывает до «пола» — погрешности сетки; дальше уменьшать шаг бесполезно.
Согласованное измельчение. Неявная схема при $\Delta t \propto \Delta x$ в итоге сходится с первым порядком: погрешность по времени убывает медленнее. Чтобы сохранить второй порядок, нужен шаг $\Delta t \propto \Delta x^2$ (как у явной схемы — но уже по точности, а не по устойчивости) или схема второго порядка по времени.

Как проверить сходимость расчёта

1
Выбрать показатели
то, ради чего строится модель: давление на забое, накопленная добыча к дате, время прорыва воды. Сходимость разных показателей различна: давления обычно сходятся быстрее, чем время прорыва.
2
Проверить шаг по времени
на рабочей сетке выполнить расчёты с шагом $\Delta t$ , $\Delta t/2$ и $\Delta t/4$ . Разности результатов показывают порядок и величину погрешности шага; при первом порядке погрешность расчёта с шагом $\Delta t/2$ примерно равна разности его результата и результата с шагом $\Delta t$ .
3
Проверить сетку
три сетки с двукратным измельчением — с шагом, погрешность которого уже мала, или с шагом, уменьшаемым вместе с ячейками. Наблюдаемый порядок по трём сеткам и экстраполяция — по формуле урока Сетки в гидродинамических моделях.
4
Сравнить порядок с теорией
наблюдаемый порядок, заметно отличающийся от теоретического, означает, что сетки ещё грубы (не в асимптотической области), решение негладкое или в расчёте есть ошибка — например, в задании проводимостей или скважин.
5
Сравнить погрешность с задачей
оценённую погрешность сравнивают с точностью, нужной для решения: с разницей между вариантами разработки, с неопределённостью исходных данных.

Расчётный примерПроверочная задача: сходимость по сетке и по шагу

Исходные данные (условные)

Закрытая полоса: длина $L$ = 2000 м; порода урока Режимы фильтрации, $\eta$ = 8412 м²/сут
Начальное возмущение: $p - p_{\text{ср}} = A_0\cos(\pi x/L)$ , $A_0$ = 1 МПа: левый конец на 1 МПа выше среднего, правый — ниже
Сравниваем: амплитуду возмущения через 30 сут на сетках из 4, 8 и 16 ячеек

Решение

Точное решение: возмущение сохраняет форму и затухает со скоростью $\lambda = \eta\pi^2/L^2 = 0{,}020756$ сут⁻¹; через 30 сут амплитуда $\e^{-0{,}6227} = 0{,}53650$ МПа. Это самая медленная составляющая выравнивания давления в закрытой полосе: любое возмущение в конце концов затухает по этому закону.
Уравнения ячеек с непроницаемыми концами сохраняют косинусную форму точно, меняется только скорость затухания: $\lambda_N = (4\eta/\Delta x^2)\sin^2(\pi/(2N))$ = 0,019711, 0,020491 и 0,020690 сут⁻¹ для $N$ = 4, 8 и 16. Погрешность сетки 5,04, 1,28 и 0,32 % — уменьшается вчетверо при каждом измельчении.
Неявная схема, шаг 3 сут (10 шагов): за шаг амплитуда умножается на $1/(1 + 3\lambda_N)$ , и через 30 сут $A_4$ = 0,56298, $A_8$ = 0,55069, $A_{16}$ = 0,54761 МПа.
Изменения при измельчении 0,01229 и 0,00308 МПа, наблюдаемый порядок $\log_2(0{,}01229/0{,}00308) = 2{,}0$ — второй, как у второй разности. Экстраполяция при $p$ = 2: $A_{16} + (A_{16} - A_8)/3 = 0{,}54658$ МПа.
Но точная амплитуда 0,53650 МПа: экстраполированное значение выше на 0,0101 МПа (1,9 %). Это погрешность шага 3 сут — сетки сходятся к решению с этим шагом.
Те же три сетки с шагом 1,5 сут: 0,55836, 0,54582, 0,54267 МПа; экстраполяция 0,54162 МПа. Погрешность шага у неявной схемы первого порядка, поэтому по двум шагам $2\cdot0{,}54162 - 0{,}54658 = 0{,}53666$ МПа — отличие от точного уже 0,00016 МПа.
Схема Кранка — Николсон с шагом 3 сут на тех же сетках: 0,55349, 0,54068, 0,53746 МПа; экстраполяция 0,53639 МПа — в пределах 0,0001 МПа от точного: погрешность её шага почти в сто раз меньше.

Ответ: Три сетки при шаге 3 сут показывают второй порядок и дают 0,54658 МПа вместо точных 0,53650 МПа; остающиеся 0,0101 МПа — погрешность шага. Учёт шага (второй ряд расчётов с 1,5 сут) приближает оценку к точной до 0,00016 МПа.

Что это значит. Проверка «только по сетке» здесь выглядела бы безупречно: порядок ровно второй, изменения убывают, экстраполяция устойчива. И всё же ответ ошибочен почти на 2 % — ровно на погрешность шага, которую сетка не видит. В реальной модели точного решения нет, поэтому шаг проверяют отдельно, как в шестом пункте. Задачи с известным решением, как эта, служат для проверки самих схем и программ — их называют верификацией, в отличие от сравнения модели с промысловыми данными (Адаптация модели к истории разработки). Данные примера условные.

РасчётПогрешность амплитуды в проверочной задаче (через 30 сут) в зависимости от шага по времени. Без погрешности сетки неявная схема сходится с первым порядком (наклон 1), схема Кранка — Николсон — со вторым (наклон 2). На сетке из 8 ячеек уменьшение шага упирается в погрешность сетки — 4,3 кПа.

Термины и обозначения

Метод конечных разностей: замена производных разностями значений в узлах сетки; превращает дифференциальное уравнение в систему алгебраических — разностную схему.
Шаблон: набор узлов сетки (по пространству и времени), входящих в разностную формулу одного уравнения.
Погрешность аппроксимации $\psi$: невязка, которую даёт точное решение, подставленное в разностное уравнение; её главный член пропорционален степени шага.
Порядок аппроксимации: показатель степени шага в главном члене погрешности; при уменьшении шага вдвое погрешность убывает примерно в $2^m$ раз.
Устойчивость: схема не усиливает ошибки от шага к шагу; вместе с аппроксимацией обеспечивает сходимость (теорема эквивалентности Лакса).
Множитель перехода $G$: во сколько раз схема меняет амплитуду гармоники за шаг; устойчивость — $|G| \le 1$ для всех гармоник.
Сеточное число Фурье $r$: $\eta\Delta t/\Delta x^2$ ; явная схема уравнения давления устойчива при $r \le 1/2$ .
Схема с весами: правая часть — взвешенное среднее значений начала и конца шага с весом конца $\sigma$ : 0 — явная, 1 — неявная, ½ — Кранка — Николсон.
Схема Кранка — Николсон: схема с весом ½: второй порядок по времени, устойчива при любом шаге, но при крупном шаге плохо гасит быстрые составляющие.

Частые ошибки и заблуждения

Неверно

«Схема второго порядка всегда точнее схемы первого».

На самом деле

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

Неверно

«Первые десятки шагов выглядят разумно — значит, схема устойчива».

На самом деле

Неустойчивая гармоника растёт из ошибок округления: при множителе −2,2 за 10 шагов в три тысячи раз, за 30 — в десятки миллиардов. Плавная часть решения при этом может считаться верно до самого «взрыва». Шаг проверяют по условию устойчивости.

Неверно

«Результат перестал меняться при измельчении сетки — значит, он точен».

На самом деле

При неизменном шаге сетки сходятся к решению с этим шагом. В проверочной задаче урока экстраполяция по трём сеткам дала второй порядок и всё равно ошиблась на 0,010 МПа — на погрешность шага 3 сут. Шаг проверяют отдельно.

Неверно

«Итерации на всех шагах сошлись — модель сошлась».

На самом деле

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

Неверно

«Чем меньше шаг численного дифференцирования, тем точнее производная».

На самом деле

Только для точных значений. При погрешности значений $\varepsilon$ к ошибке формулы добавляется ошибка порядка $\varepsilon/\Delta x$ , и при очень малом шаге шум преобладает — поэтому производную замеренного давления считают по разнесённым точкам.

Подробно: для продвинутых

Вывод погрешности аппроксимации схемы с весами

Подставим точное решение в схему и разложим все значения в ряд Тейлора около точки $(x_i,\ t^{n+1/2})$ — середины шага. Разность по времени — центральная относительно середины шага, поэтому $(p^{n+1} - p^{n})/\Delta t = \partial_t p + \tfrac{\Delta t^2}{24}\,\partial_t^3 p + \ldots$ Вторая разность в каждый момент равна $\partial_x^2 p + \tfrac{\Delta x^2}{12}\,\partial_x^4 p + \ldots$ , а значения в начале и в конце шага отличаются от значений в середине на $\mp\tfrac{\Delta t}{2}\,\partial_t$ плюс члены порядка $\Delta t^2$ .

Разложение взвешенной правой части

\sigma\,\delta^2 p^{n+1} + (1 - \sigma)\,\delta^2 p^{n} = \partial_x^2 p + \left(\sigma - \tfrac12\right)\Delta t\,\partial_t\partial_x^2 p + \frac{\Delta x^2}{12}\,\partial_x^4 p + O\!\left(\Delta t^2 + \Delta x^2\Delta t + \Delta x^4\right)

где:

$\delta^2 p$ трёхточечная вторая разность, делённая на $\Delta x^2$ , МПа/м²
$\partial_t,\ \partial_x$ частные производные по времени и координате в середине шага, —

Единицы: любые согласованные

Справедливо, если

гладкое решение, постоянные коэффициенты, равномерная сетка.

Умножим на $\eta$ и вычтем из разности по времени. Так как точное решение удовлетворяет $\partial_t p = \eta\,\partial_x^2 p$ , главные члены сокращаются, а $\eta\,\partial_t\partial_x^2 p = \partial_t^2 p$ . Остаётся $\psi = (\tfrac12 - \sigma)\Delta t\,\partial_t^2 p - \tfrac{\eta\Delta x^2}{12}\partial_x^4 p$ — формула основного раздела.

Неравномерная сетка

Пусть расстояния от узла $i$ до правого и левого соседа равны $h_+$ и $h_-$ . Обобщение трёхточечной формулы и её погрешность снова дают разложения Тейлора:

Вторая разность на неравномерной сетке узлов

\frac{2}{h_+ + h_-}\left(\frac{p_{i+1} - p_i}{h_+} - \frac{p_i - p_{i-1}}{h_-}\right) = p''_i + \frac{h_+ - h_-}{3}\,p'''_i + O(h^2)

где:

$h_+,\ h_-$ расстояния до правого и левого соседних узлов, м
$h$ характерный шаг сетки, м

Единицы: любые согласованные

Справедливо, если

сетка узлов (значения в точках); для блочно-центрированной сетки анализ сложнее, но практический вывод тот же;
гладкая функция.

Если соседние шаги различаются в фиксированное число раз, главный член пропорционален $h$ — порядок падает до первого. Если шаг меняется плавно, $h_+ - h_-$ имеет порядок $h^2$ , и второй порядок сохраняется. Отсюда правило: размеры соседних ячеек меняют постепенно, а не скачком (локальное измельчение).

Схема, которая сходится не к тому уравнению

Аппроксимацию нельзя считать очевидной. Классический пример — схема Дюфорта — Франкела (DuFort–Frankel). Возьмём явную схему с центральной разностью по времени через шаг, $(p_i^{n+1} - p_i^{n-1})/(2\Delta t)$ : для уравнения давления она неустойчива при любом шаге. Если во второй разности заменить $2p_i^n$ на $p_i^{n+1} + p_i^{n-1}$ , схема остаётся явной и становится устойчивой при любом шаге. Но разложение Тейлора добавляет к её погрешности член $\eta\,(\Delta t/\Delta x)^2\,\partial_t^2 p$ : он исчезает, только если шаг по времени уменьшают быстрее шага сетки. Если же измельчать их пропорционально, $\Delta t/\Delta x = \beta$ , схема сходится к решению уравнения $\partial_t p + \eta\beta^2\partial_t^2 p = \eta\,\partial_x^2 p$ — с «инерционным» членом, которого в фильтрации нет. Вывод: аппроксимацию проверяют для того способа измельчения, который действительно применяется.

Двумерная сетка: пятиточечный шаблон и нумерация ячеек

В двух измерениях вторые разности по $x$ и по $y$ складываются, и в уравнение ячейки $(i, j)$ входят она сама и четыре соседа по граням — пятиточечный шаблон (five-point stencil). Для явной схемы анализ гармоник даёт $G = 1 - 4r_x\sin^2(\xi/2) - 4r_y\sin^2(\zeta/2)$ с фазами $\xi$ и $\zeta$ по двум осям; худший случай — шахматная гармоника по обеим осям, откуда $r_x + r_y \le \tfrac12$ .

Для решения неявной системы ячейки нумеруют одним индексом. При естественной нумерации по строкам ячейка $(i, j)$ получает номер $m = (j - 1)N_x + i$ , где $N_x$ — число ячеек вдоль $x$ . Соседи по $x$ имеют номера $m \pm 1$ , соседи по $y$ — $m \pm N_x$ . Матрица системы содержит пять ненулевых диагоналей: главную, две соседние и две, отстоящие на $N_x$ . Ширину этой «ленты» определяет число ячеек в строке, поэтому для прямых методов исключения ячейки выгодно нумеровать сначала вдоль короткой стороны модели: у модели 150 × 30 ячеек лента сужается со 150 до 30, а объём работы ленточного исключения, пропорциональный квадрату ширины ленты, — примерно в 25 раз. Как решают такие системы с миллионами неизвестных — в уроке Решение систем линейных уравнений в симуляторах.

Связь с другими темами

Уравнение давления для однофазной фильтрации — уравнения ячеек, явная и неявная схемы, предельный шаг и метод прогонки.
Что такое гидродинамическая модель и зачем она нужна — явный и неявный шаг для одного «бака» — простейший случай анализа гармоник.
Сетки в гидродинамических моделях — проверка сходимости по сетке и экстраполяция к бесконечно мелкой сетке; численная диффузия схемы против потока.
Уравнение пьезопроводности — уравнение пьезопроводности — модельное уравнение для разностных схем.
Уравнения двухфазной фильтрации — уравнения давления и насыщенности, к которым применяют те же схемы.
Схемы IMPES и полностью неявная — схемы для многофазных нелинейных уравнений, итерации Ньютона и выбор шага.
Решение систем линейных уравнений в симуляторах — как решают разреженные системы, которые дают неявные схемы.
Численная дисперсия — погрешность схем переноса насыщенности и размытие фронтов.
Диагностика по производной давления — численное дифференцирование замеров с шумом.

Проверьте себя

1
Центральная разность при шаге 40 м ошибается на 0,6 %. Какой ошибки ожидать при шаге 10 м (функция гладкая, значения точные)?
1. А.0,6 %
2. Б.0,15 %
3. В.0,04 %
4. Г.0,01 %
Показать ответ
Правильный ответ: В. 0,04 %
У центральной разности второй порядок: погрешность пропорциональна квадрату шага. Шаг уменьшился в 4 раза, погрешность — в 16 раз: $0{,}6/16 \approx 0{,}04$ %. Вариант 0,15 % соответствовал бы первому порядку.
2
Схема с весом $\sigma$ = 0,25 решает уравнение давления на сетке с ячейками 100 м при пьезопроводности 20 000 м²/сут. При каком наибольшем шаге она устойчива и во сколько раз этот шаг больше, чем у явной схемы?
Показать ответ
Условие для $0 < \sigma < \tfrac12$ : $r \le 1/(2(1 - 2\sigma)) = 1/(2\cdot0{,}5) = 1$ , то есть $\Delta t \le \Delta x^2/\eta = 10\,000/20\,000 = 0{,}5$ сут. У явной схемы $r \le 0{,}5$ : $\Delta t \le 0{,}25$ сут. Предел вдвое больше — но и работа на шаге больше: нужно решать систему уравнений, как в неявной схеме.
3
Какая схема при крупном шаге после пуска скважины даёт затухающие колебания давления в соседних со скважиной ячейках?
1. А.Неявная
2. Б.Кранка — Николсон
3. В.Явная при $r < 1/2$
4. Г.Никакая: все устойчивые схемы монотонны
Показать ответ
Правильный ответ: Б. Кранка — Николсон
У схемы Кранка — Николсон множитель перехода самой короткой гармоники $(1 - 2r)/(1 + 2r)$ при $r \gg 1$ близок к −1: она почти не затухает и меняет знак на каждом шаге. Пуск скважины создаёт резкое изменение, богатое короткими гармониками. У неявной схемы множитель $1/(1 + 4r)$ положителен и мал. Явная схема при $r < 1/2$ к крупному шагу не относится; новое давление ячейки у неё — взвешенное среднее давлений самой ячейки и соседей с неотрицательными весами (Уравнение давления для однофазной фильтрации), поэтому давление меняется монотонно.
4
Давление на забое рассчитали на трёх сетках с двукратным измельчением при одном и том же шаге 5 сут: 14,62; 14,95 и 15,03 МПа. Каков наблюдаемый порядок, экстраполированное значение и что ещё нужно проверить?
Показать ответ
Изменения 0,33 и 0,08 МПа, порядок $\log_2(0{,}33/0{,}08) = \log_2 4{,}1 \approx 2{,}0$ — как у второй разности. Экстраполяция: $15{,}03 + 0{,}08/3 \approx 15{,}06$ МПа. Но это значение для шага 5 сут: погрешность шага сетки не видят. Нужны расчёты с шагом 2,5 сут (и 1,25 сут) на одной из сеток, чтобы оценить её.
5
Схема аппроксимирует линейное уравнение давления с погрешностью $O(\Delta t + \Delta x^2)$ и устойчива. Что можно сказать о её решении при измельчении?
1. А.Оно сходится к точному с первым порядком по времени и вторым по пространству
2. Б.Оно сходится только при $r \le 1/2$
3. В.Оно сходится, но порядок предсказать нельзя
4. Г.Оно может не сходиться: нужна ещё проверка
Показать ответ
Правильный ответ: А. Оно сходится к точному с первым порядком по времени и вторым по пространству
По теореме эквивалентности Лакса для линейной задачи аппроксимация и устойчивость дают сходимость, а ошибка решения ограничена временем расчёта, умноженным на наибольшую погрешность аппроксимации (с точностью до постоянного множителя), — значит, убывает с теми же порядками. Условие $r \le 1/2$ относится только к явной схеме; здесь устойчивость схемы дана в условии.
6
Почему исследования сходимости только по сетке при неизменном шаге по времени может быть недостаточно, даже если наблюдаемый порядок «правильный»?
Показать ответ
Погрешность складывается из пространственной $C_x\Delta x^p$ и временной $C_t\Delta t^q$ . При неизменном шаге временная часть постоянна: в разностях результатов на разных сетках она сокращается, поэтому порядок получается правильным, а экстраполяция сходится к решению с этим шагом. В проверочной задаче урока так осталась ошибка 0,010 МПа — почти 2 % амплитуды. Шаг проверяют отдельным рядом расчётов.

Обновлено: 10 октября 2026 г.