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

2025-01-02 в 8:15, admin, рубрики: дельта функция, дифференциальные уравнения, добыча нефти, преобразование фурье, свертка, уравнение теплопроводности

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

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

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

$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}$ - её производная по времени,
- коэффициент пьезопроводности,
- оператор Лапласа (мы рассмотрим для 1-мерного и 2-случаев, т.е. n=1,2),,
- внешние источники.

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

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

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

где

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

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

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

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

Свёртка

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

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

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

- единичная функция является нейтральным элементом относительно свёртки,
- тождественная функция является оператором первообразной первого порядка (по сути, это теорема Ньютона-Лейбница),
$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{2pi}} int_{-infty}^{infty} f(x)e^{-iomega x}dx \ F^{-1}[hat{f}](x)=frac{1}{sqrt{2pi}} int_{-infty}^{infty} hat{f}(omega)e^{iomega x}domega$

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

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

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

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

и, что более для нас важно, .
$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{2pi}}$

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

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

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

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

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

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

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

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

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

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

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

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

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

$u_t=frac{1}{chi} (u_{xx} +u_{yy}) + qdelta, \ -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{2pi}}.$

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

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

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

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

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

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

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

Заключение

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