Обратная задача нестационарного течения несжимаемой жидкости в трубе с проницаемой стенкой

Гидродинамические исследования движения однофазных жидкостей в перфорированных трубах помимо теоретического значения, имеют и многочисленные практические применения в пневмо- и гидротранспорте, химической технологии, разработке нефтяных и газовых месторождений горизонтальными скважинами, ракетной технике, мелиорации, гидрологии, биомеханике и т. д. [1-4]. При исследовании таких задач дискретное распределение точек перфорации (точек отбора или подкачки жидкости) заменяется непрерывным и исследования задач сводятся к изучению нестационарного движения жидкости в трубах с проницаемыми стенками. Численному и аналитическому исследованию математических моделей нестационарного течения жидкостей с различными реологическими свойствами в трубах с проницаемыми стенками посвящены работы [1—7].

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

Постановка задачи и метод решения

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

Гамзаев Х.М.                         Обратная задача нестационарного течения несжимаемой жидкости в трубе с проницаемой стенкой ^uixitl+„ (x,t) ^(У! =-1 fet) + v ^ - AuM „ (x, ,),(1)

d t            dx      p dx      dx2

и дифференциальное уравнение неразрывности потока с учетом расхода жидкости через стенки трубы:

d u ( x , t ) 4

—7—=-q(x, t)> оxd где p(x, t) - давление жидкости в трубе, u (x, t) - средняя скорость по сечению трубы, q(x, t) - скорость оттока жидкости через стенку трубы, v - кинематическая вязкость жидкости, p -плотность жидкости, X - коэффициент гидравлического сопротивления.

Предположим, что отток жидкости через проницаемую стенку трубы происходит в соответствии с законом Старлинга [5, 6], согласно которому скорость оттока жидкости прямо пропорциональна разнице между давлением жидкости внутри трубы и давлением в окружающей среде:

q (x, t) = k (t)(p (x, t) - Pe (x, t)),(3)

где p e ( x , t ) - давление внешней среды, k ( t ) - коэффициент проницаемости стенки трубы. Учитывая соотношение (3), уравнение (2) запишем в виде:

-

du (x, t) 4,... . .,

----= - k ( t )( P ( x , t ) - P e ( x , t )) .

оxd

Преобразуем систему уравнений (1), (4). Найдем p ( x , t ) из уравнения (4) , _x . d d u ( x , t ) p ^{( x} , t ) = p e ^{( x} , t )^-           ----,

4 k (t) dx и полученное выражение подставим в уравнение (1). В результате получим:

d u ( x , t )

—+ u ( x , t ) d t

d u ( x , t )      1 d p e ( x , t )

------=----+

d x

p  d x

d d ² u ( x , t )

+ v

4 p k ( t ) d x ²

d ² u ( x , t ) X u ( x , t )|

d x ²

2 d

u ( x , t ),    (6)

0 <  x <  l , 0 <  t <  T .

Предположим, что для уравнения (6) известно следующее начальное u ( x , 0) = ф ( x )

и граничные условия u (0, t) = 5( t),                                                (8)

u ( l , t ) = w ( t ).                                                    (9)

Очевидно, что, задавая функции ф(x), 5(t), w(t), k(t), решив прямую задачу (6)-(9), можно найти распределение скорости течения в трубе u (x, t). А по формуле (5) - распределение давления по длине трубы. Теперь предположим, что коэффициент проницаемости стенки трубы k(t) неизвестен и подлежит определению наряду с функцией u (x, t). Однако для корректной постановки задачи необходимо задавать дополнительное условие. Пусть в конце трубы x = l наряду со скоростью течения одновременно задается и давление жидкости p (l, t) = Y( t).

Тогда, записав уравнение (4), при x = l с учетом краевого условия для давления, в качестве дополнительного условия для уравнения (6) будем иметь du(l,t) 4 7                   w/.Ax

7----= 77k (t )(Y( t) — pe (l, t)).(10)

oxd

Таким образом, задача заключается в определении пары функций { u ( x , t ), k ( t )}, удовлетворяющих уравнению (6) и условиям (7)-(10). Поставленная задача (6)-(10) относится к классу коэффициентных обратных задач математической физики [8, 9].

Для численного решения коэффициентной обратной задачи (6)-(10), сначала построим ее дискретный аналог. С этой целью введем равномерную разностную сетку to = {(tj,xi):xi = iAx, tj = jAt, i = 0,1,2,...n, j = 0,1,2,...m},

Математика

в прямоугольной области {0 < x < l, 0 < t < T} с шагами Ax = l/n по переменной x и At = T/m по времени t . Используя полунеявную аппроксимацию по времени для нелинейных членов уравнения (6), дискретный аналог задачи (6)-(10) на сетке й представим в виде jjyj 1 jbjb jb     jib                     ?/j 1jib 1 ?/j 1 ^u,-           j              j ui   ui   + uj-1 ui   ui-1 = v ui+1   2ui + ui-1 + Jd  ui+1   2ui   + ui-1 _ I 1   I uj _ 1 pei+1   pei-1

At      i     Ax            Ax2        kj 4p       Ax2           2d   i p    2Ax i = 1,2,..., n -1, u0 = ф, i = 0, 2,..., n, u 00 = 9j, un = w , un - un-1

A x

= -k(уj -pen), d где ui = u (x- ,tj ), pei = Pe (x,tj ), kj = k(tj ), Yj = y(tj ), wj = w(tj ), 9j = 9(tj ), ф = ф(x).

Полученная разностная задача представляет собой систему линейных алгебраических уравнений, в которой в качестве неизвестных выступают приближенные значения искомых функций u ( x , t ) и k ( t ) в узлах разностной сетки < у , т. е. u- , k j , i = 0, 1, 2,..., n , j = 1, 2, 3,...., m . Данную систему разностных уравнений преобразуем к виду:

auL1-cui+ biui+1 = у-gi-1 + ft-1, i = 1,2,...,n- 1, kj u j= 9j,(12)

unj = wj ,(13)

un = un-1 -kj 4(Yj -Pejn)Ax,(14)

nd j = 1,2,..., m, u0 = ф,, i = 0, 2,..., n,(15)

где

^ | u,j 1 | a, = vAt + uj AtAx, b = vAt, c,- = a,- + b + Ax +--i Ax At, i         i          ii ii

2d j-1     j-1j-Цл/          tn

„j - 1 _ d ⁽ u i + 1 ² ui ⁺ u i - 1 ^{) A} t r j - 1 _ ⁽ p ei + 1 p ei - 1^{) A} t ^A x a_v2, . - - 1

^gi =               .                , Ji = _n               ^{A x u}i .

4 P                    2 p

С целью разделения разностной задачи (11)–(15) на взаимно независимые подзадачи, каждая из которых может решаться самостоятельно, решение этой системы при каждом фиксированном значении j , j = 1,2,..., m представим в виде [10, 11]:

uj = ^ +—nj, i = 0,1,...,n, kj где ^i, ni — неизвестные переменные. Подставляя выражение для uj в (11)-(13), получим

Г а й - 1 - е й + b ^ + 1 - ff¹ 1 + J- [ « i S - - c^i + Ьйи - g i ^{- 1} 1 = 0, S ' + ^A S = 9 , S + J - n ', = w .

L                                             —I I^J L                                            —I                  I^JI^J

Из последних соотношений получим следующие разностные задачи для определения вспомогательных переменных Si , П aiSi- -cSj + bS+1 = Jij-1,(17)

Sj= 9j ,(18)

Sn = wj .(19)

ain/-1- cni+bin/+1=gi-1,

П 0 ^— О,

ni ^— 0.

Разностные задачи (17)-(19) и (20)-(22) при каждом фиксированном значении j = 1, 2,..., m представляют собой систему линейных алгебраических уравнений с трехдиагональной матрицей, и решения этих систем можно найти методом Томаса [9].

Подставив представление (16) в (14), будем иметь:

1                     14

& + ту n nj = SU + 77 п 3 — 1 ^- k ^j - ⁽Y ^{j -} P en )^ x .

kj             kj

Отсюда можно определить приближенное значение коэффициента проницаемости стенки трубы k ^j как положительный корень данного уравнения

, j _ ^- d ⁽Y ^{- j} ) ± V d ²⁽ ^ ^{- j )2} - ¹⁶ d A x ( П 3 - п 3 — 1 ⁾⁽ Y ^{j -} Pin )

k                                                                            .(23)

4( Y - P ^jn ) A x/d

Определив k ^j по формуле (23), можно последовательно найти u ( , u ^j ,...,u³_{n - 1} по рекуррентной формуле (16).

После определения распределения скорости жидкости по длине трубопровода можно перейти к определению распределения давления. Построив дискретный аналог уравнения (5) на сетке to, получим следующую расчетную формулу для вычисления давления pi — pJjt- .d,- ui+1Л ui ,1 — 012,...,3- 1.                          (24)

• 4 k ^j   A x

Таким образом, вычислительный алгоритм решения коэффициентной обратной задачи (6)-(10) по определению коэффициента проницаемости стенки трубы k(t) при каждом дискрет ном значении временной переменной tj, j —1,2,..., m основан на:

• 1)    решении двух линейных разностных задач второго порядка (17)-(19) и (20)-(22) относительно вспомогательных переменных Y , n i , i — 1, п ;

• 2)    определении к³ из(23);

• 3)    использовании представления (16) для вычисления u i , i — 1, п - 1.

Результаты численных расчетов

Предложенный вычислительный алгоритм был опробован на модельных задачах. Численные расчеты проводились по следующей схеме:

- для заданных функций ф (x ), 0 (t ), w ( t ), к ( t ) определяется решение прямой задачи (6)-(9), т. е. функция u ( x , t ), 0 <  x <  l , 0 <  t <  T ;

- найденная зависимость

d д u ( l , t )

p(l, t) — p (l, t)----—- — Yt) принимается в качестве точ- e       4 к (t) дx ных данных для решения обратной задачи по восстановлению к(t).

Численные расчеты выполнялись на разностной сетке с шагами A x = 0,5 м, A t = 0,01 с. Результаты численного эксперимента, проведенного для случая d = 0,4 м, l = 100 м, У ( t ) = 1,5 м/с, w ( t ) = 0,5 м/с, ф ( x ) — 1,5 - x11 м/с, p = 1000 кг/м³, v = 10^-6 м²/с, к ( t ) = 0,05 м²-с/кг, к ( t ) — 0,4 - 0,3sin10 1 м²-с/кг, к ( t ) — 0,1/V t м²-с/кг, X = 0,02 представлены в таблице.

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

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