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

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

Уравнение пьезопроводности выглядит следующим образом:

$\frac{\partial u}{\partial t} = \frac{1}{\chi} \Delta u + f(x,t), \\ u(x,0)=0, \quad x \in R^n, t > 0,$

где

- скачок давления (наша искомая функция),
$\frac{\partial u}{\partial t}$ - её производная по времени,
$\chi > 0$ - коэффициент пьезопроводности,
$\Delta$ - оператор Лапласа (мы рассмотрим для 1-мерного и 2-случаев, т.е. n=1,2),,
- внешние источники.

Отметим, что слово "скачок" в определении функции очень важно. Дело в том, что нефть в пласте перед началом добычи находится под давлением , которое во много раз может превышать атмосферное. Поэтому удобно рассматривать модель с точки зрения именно скачка относительно начального давления. Исходное давление связано с этим скачком по формуле

Теперь перейдем ко внешним источникам. Естественно, в нашей постановке внешним источником является скважина. Моделировать её можно по-разному, например можно считать

$f(x,t) = \cases{q, \text{ если } ||x - x_w|| <= r_w, \\ 0, \text{ иначе};}$

где

$q \in R$ - дебит скважины - некоторая константа, которая определяет силу, с которой закачивается или выкачивается жидкость,
- координата центра скважины,
- радиус скважины.

Мы же рассмотрим предельный случай этой модели при $r_w \rightarrow 0$ . Именно здесь и возникает особенность задачи, которую называют точечным источником. Математически записать это можно с помощью дельта-функции Дирака:

$f(x,t) = q \delta(x-x_w).$

Для простоты будем считать, что скважина находится в начале координат, т.е. а производные будем записывать в индексе. Тогда для 1-мерного и 2-мерного случаев получим:

$u_t = \frac{1}{\chi}u_{xx} + q \delta, \\ u_t=\frac{1}{\chi}(u_{xx}+u_{yy})+q\delta.$

Свёртка

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

$f \ast g= \frac{d}{dt} \int_0^t f(\tau)g(t-\tau)d\tau.$

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

$f \ast 1 = f$ - единичная функция является нейтральным элементом относительно свёртки,
$f' \ast t = f(t) - f(0)$ - тождественная функция является оператором первообразной первого порядка (по сути, это теорема Ньютона-Лейбница),
$t \ast ... \ast t \text{ (n раз) } = \frac{t^n}{n!}$ - первообразная n-го порядка это свёртка тождественных функций.

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

$e^{-at} \ast t = \frac{e^{-at}-1}{-a} \Longrightarrow ae^{-at} \ast t = 1-e^{-at},$

Операция свёртки линейна, поэтому мы можем собрать экспоненты в левой части равенства и вынести за скобки:

$e^{-at} \ast (a t + 1) = 1,$

Мы получили не что иное, как функцию $f^{-1}(t)=at+1$ обратную к экспоненциальной функции $f(t)=e^{-at}$ , для которых выполняется равенство:

$f \ast f^{-1}=1.$

Давайте для примера решим ДУ

$y'+ay=g(t),\\ y(0)=y_0.$

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

$y'+ay=g(t) \quad | \ast t \quad \\ y-y(0) + a y \ast t = g \ast t \\ y \ast (1 + at) = y_0 + g \ast t \\ y \ast (1 + at) = y_0 + g \ast t \quad | \ast e^{-at} \\ y = (y_0 + g \ast t) \ast e^{-at}.$

И вот мы уже имеем ответ $y(t) = (y_0 + g \ast t) \ast e^{-at}$ . К сожалению, такой способ не всесилен, поскольку решение можно получить лишь тогда, когда коэффициенты уравнения постоянны, а решить с переменными возможно только при особых случаях.

Преобразование Фурье

Для того, чтобы получить аналитическое решение нам также потребуется преобразование Фурье. Как известно, одномерное прямое и обратное преобразования определяются по формулам

$\hat{f}(\omega)=F[f](\omega) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} f(x)e^{-i\omega x}dx \\ F^{-1}[\hat{f}](x) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} \hat{f}(\omega)e^{i\omega x}d\omega$

Функции и $\hat{f}$ при этом называются оригиналом и изображением соответственно.

Также нам потребуются преобразования в двумерном случае:

$F[f](\omega_1, \omega_2) = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}^2} f(x)e^{-i(\omega_1 x + \omega_2 y)}dxdy \\ F^{-1}[\hat{f}](x, y) = \frac{1}{\sqrt{2\pi}} \int_{\mathbb{R}^2} \hat{f}(\omega_1, \omega_2)e^{i(\omega_1 x + \omega_2 y)}d\omega_1 d\omega_2$

В общем-то, эти формулы нам не нужны, достаточно его свойств и пары равенств из таблицы преобразований.

$F[f'] = i \omega F[f]$ и, что более для нас важно, $F[f''] = - \omega^2 F[f]$ .
$F[e^{-ax^2}](\omega)= \frac{1}{\sqrt{2a}} e^{-\omega^2 /(4a)}$ или $F^{-1}[e^{-\omega^2/(4a)}](x)= \sqrt{2a} \ e^{-ax^2}$ .
$F[\delta](\omega) = \frac{1}{\sqrt{2\pi}}$

Решение при n=1

Уравнение будет иметь вид

$u_t = \frac{1}{\chi} u_{xx} + q\delta, \\ -\infty < x < +\infty, t > 0, \\ u(x,0)=0.$

Применим преобразование Фурье:

$\hat{u}_t + \frac{\omega^2}{\chi} \hat{u} = \frac{q}{\sqrt{2\pi}}.$

Мы получили обыкновенное дифференциальное уравнение относительно переменной и изображения $\hat{u}$ . Решим его с помощью свертки:

$\hat{u} = \frac{q}{\sqrt{2\pi}} t \ast e^{-\omega^2 t/\chi} .$

И сейчас важный момент! Мы могли бы вычислить свёртку и потом получить ответ, применив обратное преобразование Фурье. Но тогда мы получили бы функцию, преобразования которой нет в таблице. Поэтому мы сделаем наоборот: сначала обратное преобразование Фурье, и потом свёртка. Мы имеем право так сделать, поскольку свёртка идет по переменной , обратное преобразование - по $\omega$ , и, к тому же, $\omega$ встречается лишь во втором операнде. Ещё одно удобство свёрточной записи, которое мы могли бы не заметить, если решали бы ДУ обычным методом!

Считая, что $\frac{1}{4a} = \frac{t}{\chi} \Rightarrow a = \frac{\chi}{4t}$ , применим обратное преобразование Фурье:

$u= q \sqrt{\frac{\chi}{4\pi}} \ t \ast \frac{e^{-x^2 \chi/(4t)}}{\sqrt{t}}= q \sqrt{\frac{\chi}{4\pi}} \int_0^t \frac{e^{-x^2 \chi/(4\tau)}}{\sqrt{\tau}} d\tau.$

Как хорошие математики, мы должны были бы проверить корректность полученного решения, подставив в исходное уравнение, но сделать это так просто не получится. Полученное нами выражение не является элементарной функцией. Если же действительно хочется проверить решение, то можно привести его к выражению, которое содержит, например, $\operatorname{erfc}(x)$ - дополнительную функцию ошибок. Однако вывод этого выражения очень сложный и громоздкий, а моя же цель не состоит в том, чтобы нагружать вас тяжелой математикой и длинными формулами. Поэтому мы примем это на веру и пойдем дальше.

Решение при n=2

Уравнение будет иметь вид

$u_t = \frac{1}{\chi} (u_{xx} +u_{yy}) + q\delta, \\ -\infty < x,y < +\infty, t > 0, \\ u(x, y, 0)=0.$

Аналогично, применим двумерное преобразование Фурье (можно представить как применение одномерного по переменным и ):

$\hat{u}_t + \frac{\omega_1^2 + \omega_2^2}{\chi} \hat{u} = \frac{q}{\sqrt{2\pi}}.$

Решаем при помощи свёртки:

$\hat{u} = \frac{q}{\sqrt{2\pi}} t \ast \operatorname{exp}({-\frac{\omega_1^2 + \omega_2^2 }{\chi} t}) .$

Здесь всё также $a = \frac{\chi}{4t}$ . После обратного преобразования Фурье (как и ранее - двойное одномерное по $\omega_1$ и $\omega_2$ ) получим:

$u = \frac{q}{2\pi} \ t \ast \frac{\chi}{2t} \operatorname{exp}(-\chi \frac{x^2+y^2}{4t}) = \frac{q\chi}{4\pi} \int_0^t \frac{1}{\tau} \operatorname{exp}(-\chi \frac{x^2+y^2}{4\tau})d\tau$

Данный интеграл тоже не является элементарной функцией, но её можно записать в более удобном виде:

$u = \frac{q\chi}{4\pi} \operatorname{Ei}(-\chi\frac{x^2+y^2}{4t}),$

где $\operatorname{Ei}(x) = -\int_{-x}^{\infty} \frac{e^{-x}}x dx$ — интегральная показательная функция.

Заключение

Вот мы и получили аналитические решения для уравнения пьезопроводности с точечным источником. К сожалению, эти решения получены с учетом многих допущений. В реальных же задачах не бывает нефтеносных пластов с бесконечной границей, количество точечных источников может исчисляться десятками и сотнями, да и сам дебит может изменятся во времени, а не быть постоянной. На самом деле, даже решения этих уравнений не нужны - намного важнее находить его параметры, как например, коэффициент $\chi$ . В реальности нужно определять свойства пласта, исходя из его поведения в настоящем, т.е. решать так называемую обратную задачу, чтобы прогнозировать объем добычи в будущем. Поэтому численные методы в большинстве случаев бывают полезнее аналитических. Тем не менее, такие идеализированные модели по прежнему полезны для изучения, поскольку позволяют нам лучше понять процессы, происходящие в пласте, и находить новые приёмы для их моделирования.

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