О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ТЕПЛОВОГО ВОЗДЕЙСТВИЯ НА ПРИЗАБОЙНУЮ ЗОНУ ПЛАСТА Текст научной статьи по специальности «Математика»

УДК 622.276:536.3:519.6

О ЧИСЛЕННОМ МОДЕЛИРОВАНИИ ТЕПЛОВОГО ВОЗДЕЙСТВИЯ НА ПРИЗАБОЙНУЮ ЗОНУ ПЛАСТА

Х. М. Гамзаев, С. О. Гусейнзаде, Г. Г. Гасымов

Азербайджанский государственный университет нефти и промышленности, [email protected], [email protected], [email protected]

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

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

ABOUT NUMERICAL MODELLING OF THERMAL IMPACT ON THE BOTTOMHOLE ZONE OF LAYER

Kh. M. Gamzaev, S. O. Guseinzade, G. G. Gasymov

Azerbaijan State Oil and Industry University, [email protected], [email protected], [email protected]

The article considers a process of heat propagation from a heated gallery of production wells to an oil reservoir. A one-dimensional diffusion-convection equation for a region with an unknown moving boundary is proposed to describe this process. For the correct formulation of the problem, an auxiliary condition is set for the temperature at the gallery of wells. By changing variables, the problem is transformed into a coefficient inverse problem. Further, the time derivative is discretized and explicit-implicit schemes are used to approximate the operators of the problem. To solve the resulting differential-difference problem, a special representation is proposed. As a result, at each discrete value of the time variable, the differential-difference problem splits into two direct boundary-value problems and a linear algebraic equation with respect to the approximate value of the sought-for coefficient. The Thomas algorithm is used for the numerical solution of the obtained direct boundary value problems. Numerical experiments were performed on the basis of the proposed numerical algorithm.

Keywords: thermal treatment of formation, boundary value problem with a moving boundary, rectification method of fronts, coefficient inverse problem, explicit-implicit schemes, differential-difference problem.

Введение. Известно, что при эксплуатации нефтяных скважин в призабойной зоне образуются органические и неорганические структуры, приводящие к повышению фильтрационного сопротивления и снижению производительности скважин. Для снижения фильтрационного сопротивления и восстановления потенциальной производительности скважин применяются различные методы воздействия на призабойную зону. В зависимости от механизма, обусловливающего улучшение фильтрационных свойств призабойной зоны, различают химические, физические и тепловые методы воздействия. В настоящее время важная роль отводится тепловым методам воздействия на призабойную зону. Тепловое воздействие на призабойную зону пласта позволяет значительно снизить вязкость нефти, вернуть в растворенное состояние выпавшую из нефти твердую фазу, снижает вероятность образования парафинистых и смолистых отложений, существенно влияет как на текущие фильтрационные характеристики, так и на конечную нефтеотдачу [1-3].

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

Постановка задачи и метод решения. Пусть рассматривается прямолинейно-параллельный приток пластовой нефти к галерее эксплуатационных скважин, расположенной в сечении х = 0 прямоугольного пласта постоянной толщины и ширины. Тепловое воздействие на пласт осуществляется путем электропрогрева галереи эксплуатационных скважин. Из-за теплового влияния нагреваемой галереи скважин в пласте образуется область с подвижной границей s(t), причем ,'(г) > 0, г > 0, т. е. область теплового влияния галереи в пласте монотонно расширяется. Предполагается, что кровля и подошва пласта теплоизолированы, а теплообмен между скелетом пласта и фильтрующейся сквозь него нефтью совершается мгновенно, вследствие чего их температуры равны. Учитывая, что перенос тепла в пласте осуществляется, во-первых, конвективным путем, т. е. вместе с движущимися частицами нефти и, во-вторых, за счет теплопроводности вдоль пласта, математическую модель процесса прямолинейно-параллельного распространения тепла от нагреваемой галереи скважин в пласт можно представить в виде

с^ Щ,- • (х,,) .П, = (> < х < £(,),0 <, *,'}. (1)

дг дх ох

где Т (х, г) - температура в пласте; Л - коэффициент теплопроводности пласта; с = тс^ + (1 - т)ср, с^, ср - объемные теплоемкости нефти и скелета пласта; т - коэффициент пористости; у(г) - скорость притока нефти к галерее скважин.

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

начальный момент времени г = 0 распределение температуры в пласте и положение подвижной границы известны, т. е. для уравнения (1) имеем следующие начальные условия

т(х,0) = т, (2)

5(0) = 0. (3)

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

+0 = е(£). (4)

Предполагая, что температура на подвижной границе равна начальной температуре в пласте Т5, условие на подвижной границе можно представить в виде

т (з(г), г) = Т . (5)

Необходимо отметить, что уравнение (1) выполняется в области с подвижной границей, а закон перемещения подвижной границы неизвестен. Следовательно, для корректной постановки задачи необходимо задавать дополнительное условие. Учитывая, что закон изменения во времени температуры на галерее скважин также известен, то в качестве дополнительного условия можно принимать

Т (0, г) = Т№ (г). (6)

Таким образом, задача определения температурного поля пласта при тепловом воздействии на призабойную зону состоит в определении функций Т(х, г), 5 (г), удовлетворяющих уравнению (1) с заданными коэффициентами су, ср, с,Л, у(г) и заданным условиям (2)-(6). Существенной особенностью данной задачи является наличие подвижной границы, закон перемещения которой определяется в ходе решения задачи. Задача (1)-(6) относится к классу краевых задач с подвижной границей [4-5].

Используя идею метода выпрямления фронтов, преобразуем задачу (1)-(6). Введем замены переменных

х

у = —, г = г, Т(х, г) = Т(у, г).

5(г)

Тогда область задания уравнения (1) £28 отображается на область £ = ¡0 < у < 1, 0 < г < г * } и задача (1)-(6) в переменных (у, г) принимает вид

с52(,)Щуй = +(+ <М0 ^ , (у, О-£ = ¡0 < у < 1,0 < г < г* },(7)

дг ду \ йг ) ду ^ '

Т(у,0) = Т, (8)

5(0) = 0, (9)

- Л дТ(0, г) + с гу(г>(г)Т(0, г) = е(г>(г), (10)

ду

Т(1,г) = Т, (11)

Т (0, г) = Т№(г). (12)

Как видно из (7)-(12), в результате преобразования задача (1)-(6) с подвижной границей преобразуется к обратной задаче по определению коэффициента ¿(г). Преимущество такого преобразования заключается в том, что полученная коэффициентная обратная задача (7)-(12) рассматривается в прямоугольной области О с фиксированными границами.

Отметим, что вопросы существования, единственности и разрешимости, а также некоторые подходы к численному решению коэффициентных обратных задач для параболических уравнений исследованы в [6—11].

Для численного решения полученной коэффициентной обратной задачи (7)-(12) используем подход, предложенный в [12]. Сначала введем равномерную разностную сетку в области

0 < г < г по переменной г

Щ = \[] = ] Л ] = 0, т}

л г ди йя . -—

с шагом Лг = —. Производные — и — в уравнении (7) при г^, ] = 1, т дискретизируем разно-

т дг йг у

стью «назад»

ди( х, г)

дг

и (х, г у)- и(х, г}_1) йз _ ¿(г у) - ¿(г-)

Лг ' йг 4 * Лг

Используя явно-неявные схемы для аппроксимации дифференциальных операторов в (7) и (10), задачу (7)-(12) запишем в следующем виде

фу-1)

у-1,2 Т3 (У) - Т]-1(У) г й Т (у) Г

Л

йу2

- +

+ суз

3-1

Лг

йТ]-1( у)

йу

-ЛйТ(0) + с^ЗТ1-1 (0) = е]з],

йу

Т](1) = Т, Т](0) = т],

] = 1,2,..., т,

Т 0(у) = Т, 0 < у < 1, з0 = 0 , где Т](у) * Т(у,г]), * з(£]), V = у(г]), Т] = Тк] е] = е(г }).

0 < у < 1,

(13)

(14)

(15)

(16)

(17)

Предположим, что решение задачи (13)—(17) на каждом временном слое ] = 1, 2, ... , т можно представить в виде

Т](у) = и (у) + Зр1(у), (18)

где и1 (у), р} (у), ] = 1, т — неизвестные функции. Подставив соотношение (18) в (13)—(15), бу-

дем иметь

(у)-Т]-1 (у)_я йV(у) 1 су]2 йТ]\у)

Лг

йу2

+ -

Лг

йу

+

-1)2

р] (у ,й2Р(у) суз]-1 йТ]-1(у)

]-1 АТ ]-11

+

Лг

-Я

йу2

Лг

йу

йТ]-1( у) йу

= 0,

х,г

и

. ё]ы(0) Л-— + 5

йу

+, уТ-1(0) -е

йу у

= 0,

и (1)+(1) = т.

Из последних соотношений можно получить следующие прямые краевые задачи относительно функций и1 (у), р {у)

ф}-1у и (у) - Т1-1 {у) л йУ(у) | су (5]-1) йТ-1( у) =0

Аг

йу

Аг йу

= 0,

йу

и1 (1) = Т.

0 < у < 1,

(19)

(20) (21)

ф7-1 )2

р1 (у .й2 р'(у) у1-1 йТ1 -1(у)

-1 -1

Аг

— Л

йу2

Л

Аг йр1(0)

йу

йу

- с V

йТ1 -1( у) йу

= 0,

0 < у < 1,

+ с/у1Т1-1(0) - е1 = 0,

(22)

(23)

(24)

р1 (1) = 0 .

1 = 1,2,..., т, Т0(у) = Т, 0 < у < 1, 50 = 0. Подстановка (18) в дополнительное условие (16) дает

и (0) + 5^(0) = Р. (25)

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

Т} (у), 5 , ] = 1,2,..., т :

- для фиксированного значения временного слоя ] определяются решения задач (19)—(21) и (22)—(24), т. е. функции и1 (у) и р (у) на отрезке [0,1];

- из соотношения (25) определяется приближенное значение искомой функции ф) при

г = г1

1_Тч - и1 (0)

=

р (0)

- по формуле Т} {у) = и1 (у) + 5р1 (у) определяется приближение искомой функции Т(у, г) при г = г];

- при переходе на следующий временной слой описанная процедура вычислений снова повторяется.

Для численного решения задач (19)—(21) и (22)—(24) можно использовать метод конечных разностей. Введем равномерную разностную сетку в области [0 < у < 1] по переменной у

а у = [уг = гАу, г = 0, п, Ау = 1 / п).

Дискретные аналоги задач (19)—(21) и (22)—(24) на сетке а представим в виде

2

] , 2 и] - Т/ 1 и]+1 - 2и/ + и- , су(з]1) Т+1 - Т]

с(з] )2 —----Я-Ш--^+ ^----'— = 0 , г = 1,п-1,

Лг Лу2 Лг Лу

= 0, (26) Лу ' '

ип = Т.

с(з]-1)2 р- -Я Лг

]-К2 Р- ,Р/+1 - 2Р- + Р/-1

Лу2

Ссуз]-1 , —-+ сV

-1 -1 Т-+1 - Тг

V

Лг

V

У

Лу

= 0,

Я Р Р0 + с^ТГ1 - е] = 0,

Лу

уг0 _ уг

рП = 0.

] = 1, 2,..., т,

0 < - < п , з0 = 0 .

г = 1, п -1,

(27)

где и- * и (у- X т]1 * т]-1 (у- ), рг * р] (у- ) .

Разностные задачи (26) и (27) при каждом фиксированном значении ] = 1, 2,..., т представляют собой систему линейных алгебраических уравнений с трехдиагональной матрицей. Для решения этих систем можно использовать алгоритм Томаса (метод прогонки).

Результаты численных расчетов. Для выяснения эффективности практического применения предложенного вычислительного алгоритма были проведены численные эксперименты для модельных задач. Схема численного эксперимента заключалась в следующем. Для заданных функций е(г), з(г) и Т решалась прямая задача (7)—(12). Найденная зависимость Т(0, г) принималась за точные данные для численного решения обратной задачи по восстановлению з(г) . Результаты численного эксперимента, проведенного для случая Я = 1.28 Вт/м.С, с = 1378000 Дж/м3.С, в; = 1875 Дж/м3.С, е = 1200 Дж/м2С, Т = 20 С, у(г) = 2-10-4 м/сек, з(г) = 2-10-3 г2 м представлены в таблице. Результаты численных экспериментов показывают, что искомая подвижная граница з(г ) восстанавливается с достаточно высокой точностью при всех расчетных сетках по времени. При этом максимальная относительная погрешность восстановления искомой функции не превышает 0,052 %.

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

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

Таблица

Численные результаты по задаче

Значение функции ), м

t, сек

У точное вычисленное при А = 0.1 сек, Ах = 0,05 м

10 0,2 0,2

20 0,8 0,8

30 1,8 1,8

40 3,2 3,199

50 5,0 4,999

60 7,2 7,198

70 9,8 9,797

80 12,8 12,798

90 16,2 16,197

100 20,0 19,994

110 24,2 24,189

120 28,8 28,889

130 33,8 33,798

140 39,2 38,196

150 45,0 44,984

160 51,2 51,196

170 57,8 57,769

180 64,8 64,798

190 72,2 72,196

200 80,0 79,967

Литература

1. Сучков Б. М. Температурные режимы работающих скважин и тепловые методы добычи нефти. М. : ИКИ, 2007. 406 с.

2. Басниев К. С., Власов А. М., Кочина И. Н., Максимов В. М. Подземная гидравлика. М. : Недра, 1986. 303 с.

3. Бурже Ж., Сурко П., Комбаржу М. Термические методы повышения нефтеотдачи пластов. М. : Недра, 1988. 424 с.

4. Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача. М. : Едиториал УРСС, 2003. 758 с.

5. Гамзаев Х. М. Численное решение задачи ненасыщенной фильтрации с подвижной границей // Электрон. моделирование. 2015. Т. 37. № 1. С. 15-24.

6. Иванчов Н. И., Побыривска Н. В. Об определении двух зависящих от времени коэффициентов в параболическом уравнении // Сиб. матем. журн. 2002. № 43:2. С. 406-413.

7. Камынин В. Л. Обратная задача определения младшего коэффициента в параболическом уравнении при условии интегрального наблюдения // Матем. заметки. 2013. Т. 94. Вып. 2. С. 207-2175.

8. Костин А. Б. Восстановление коэффициента перед ut в уравнении теплопроводности по условию нелокального наблюдения по времени // Журн. вычислит. математики и матем. физики. 2015. Т. 55. № 1. С. 89-104.

9. Кожанов А. И. Параболические уравнения с неизвестными коэффициентами, зависящими от времени // Журн. вычислит. математики и матем. физики. 2017. № 6. С. 961-972.

10. Yang L., Yu J.-N., Deng Z.-Ch. An inverse problem of identifying the coefficient of parabolic equation // Applied Mathematical Modelling. 2008. V. 32. Is. 10. P. 1984-1995.

11. Kerimov N. B., Ismailov M. I. An inverse coefficient problem for the heat equation in the case of nonlocal boundary conditions // Journal of Mathematical Analysis and Applications. 2012. V. 396. Is. 2. P. 546-554.

12. Гамзаев Х. М. Численный метод решения коэффициентной обратной задачи для уравнения диффузии - конвекции - реакции // Вестн. Томск. гос. ун-та. Математика и механика. 2017. № 50. С. 67-78.