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

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

Д Р = X^- l , 2 d ’ где Д Р - перепад давления на участке трубопровода длиной l , d - диаметр трубопровода, X -коэффициент гидравлического сопротивления, р - плотность жидкости, и - средняя скорость по сечению трубопровода.

Данную формулу, а также явное выражение для коэффициента гидравлического сопротивления для ламинарного режима можно получить из точного решения уравнения стационарного течения однородных несжимаемых жидкостей по трубопроводу при соответствующих реологических законах. При этом в качестве граничного условия на стенке трубопровода используется так называемое «условие прилипания». Так возникает известный параболический профиль скорости в стационарных течениях вязких жидкостей под действием перепада давления. Однако многие исследователи на основании молекулярных гипотез приходят к выводу, что вместо условия прилипания на твердой стенке трубопровода имеет место условие скольжения [7, 8]. В литературе рассматриваются три модели взаимодействия жидкостей с твердой стенкой, которые соответствуют следующим граничным условиям: прилипание, проскальзывание по закону Навье и проскальзывание с предельным напряжением [9-12]. Однако при моделировании течения жидкостей в трубопроводах практически невозможно определить, какое из этих граничных условий реализуется на стенке трубопровода. В связи с этим для практики трубопроводного транспорта важное значение имеет исследование по определению закона распределения скоростей в поперечном сечении трубопровода для нестационарных потоков транспортируемых вязких жидкостей при неизвестном условии на стенке трубопровода.

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

Пусть имеется горизонтально расположенный простой трубопровод с жесткими стенками, длиной l, радиусом R, и по нему перекачивается вязкая несжимаемая жидкость. Предполагается, что ось 0z направлена по оси трубопровода и поток направлен вдоль оси трубы так, что из трех компонент скорости (ur, иф, uz ) остается лишь одна uz ^ 0, а ur = 0 и иф = 0. Считая течение жидкости осесимметричным, полную систему дифференциальных уравнений, описывающих данное течение, можно представить в виде [4]

∂u     ∂u ν ∂   ∂u    1 ∂P z +u z = r z -      , 0

^∂uz=0, ^∂uz=0, ^1∂P=0, ^1∂P=0 ∂z ∂ϕ ρ ∂r ρ ∂ϕ

где uz - компонента скорости течения жидкости, направленная параллельно оси трубопровода, Р - давление, р - плотность жидкости, v = Цр - кинематическая вязкость, ц - динамическая вязкость жидкости. Из второго и третьего уравнения системы (1) следует, что uz представляет функцию только r и t, а из двух последних - независимость давления Р от r и ф. А это означает, что др является функцией только времени.

Полагая

∂P ∆P(t)

u(r,t) = u_z(r,t), -    =        ,

∂z      l из системы (1) придем к следующей форме уравнения нестационарного течения вязкой несжимаемой жидкости по трубопроводу

∂u ν∂ ∂u    1∆P(t)

= r +         , 0

Плотность жидкости р, кинематическая вязкость v и перепад давления ДР(t) считаются заданными.

Пусть для уравнения (2) задаются начальное условие u |t=0=ψ(r) ,                                                (3)

и естественное граничное условие ограниченности решения при r = 0, которое эквивалентно условию

∂u

|_r=0= 0.                                                (4)

∂r

Однако в связи с тем, что скорость течения жидкости на стенке трубопровода не доступна непосредственному измерению и не может регулироваться, сформулировать граничное условие, соответствующее взаимодействию жидкости с твердой стенкой трубопровода не представляется возможным. Следовательно, чтобы сформулировать корректную задачу, помимо условий (3), (4) необходимо задавать дополнительное условие. Предположим, что закон изменения во времени объемного расхода жидкости в трубопроводе Q(t) известен. Тогда дополнительное условие для уравнения (2) можно представить в виде

R j2nrudr = Q(t).                                        (5)

При этом предполагается выполнение условия согласования

R

2л|и r) dr = Q (0).

о

Таким образом, задача заключается в определении в прямоугольной области {0 < г< R, 0 < t< T} функции u(r, t), удовлетворяющей уравнению (2) и условиям (3)-(5). Задача (2)-(5) относится к классу нелокальных задач с интегральными условиями для дифференциальных уравнений в частных производных [13-16]. Необходимо отметить, что задание нелокального условия (5) приводит к возникновению значительных трудностей при численном решении нелокальной задачи (2)-(5) именно из-за отсутствия какой-либо информации об искомом решении на границе области.

Метод решения

Задачу (2)-(5) сведем к задаче с локальными условиями. Дифференцируем соотношение (5) по переменной t

_ R du , dQ

2n r—dr = —.

0 d t      dt

Подставим в это соотношение выражения для du из уравнения (2)

„ R Гц д - дu) 1 ДР(t) 1 , dQ 2п г--1 r— I +--— dr = —.

0 L r д r V д r) p l J dt

Выполнив интегрирования по частям и учитывая условие (4), получим д u . R2 ДР (t) = 1 dQ дг r=R 2p l 2n dt

Разрешив последнее уравнение относительно дЦ- r=_R , получим недостающее граничное условие на стенке трубопровода для уравнения (2)

д u       1 dQ R

— | _r=_R =-----ДР (t).                             (6)

дr     2nvR dt 2plv

Теперь задача заключается в определении функции u(r, t), удовлетворяющей уравнению (2) и начальному условию (3) и локальным граничным условиям (4), (6).

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

- в прямоугольной области {0 < r < R, 0 < t < T} с шагами Дr = R/n по переменной r и Дt = T/m по времени t. Пользуясь интегральным методом дискретизации, уравнению (2) во внутренних узлах сетки to поставим в соответствие неявную разностную схему [ 17]

j+1 ui

^^^^^^^е

^{u j}

д t

v гД r

r+1/2

j+1 ^ui+1

^^^^^^^e

j+1 u j +1

Д r

^^^^^^^e

r-1/2

j+1 u j +1

^^^^^^^e

u

. j+1

i-1

Д r

+—др^j+1, pl

где wj = u (r, tj), ДР⁷ = ДР (tj), r±1/₂= г + Д r/2, i = 1,2,3,..., n -1, j = 1,2,3,...., m -1.

Начальное условие (3) аппроксимируется точно:

ui³ = Yi, i = 0, n.

Для повышения порядка аппроксимации граничных условий (4), (6) снова используем интегральный метод дискретизации. В результате будем иметь:

^u0    ^u0 ^r1/2

Д t    2

= ^Vr1/2

uj⁺¹- u j⁺¹

Д r

+ -¹²- ДР^j+1 , 2 pl

j +1     j 2    2               j +1     j              j +1     j +1     2

un    un ^R   rn-1/2 _ ^{1 Q}Q           un    ^un-1   ^rn-1/2 A +1

—                V r—A / Д                   ДР

Д t        2       2п    Д t        n 1/2    Д r       2 pl где Yi = Иr), Q =Q(tj )•

Преобразуем полученную систему разностных уравнений к следующему виду: au+¹- cu⁺¹+ bu-ij+A = -( u +—APj⁺¹), i = 1,2,3,.... n -1, j = 1,2,3,...., m -1, I i—1       II I i+1           I

pl ___ ui = ¥i, i = 0, n, u j+1 = 61 uj+1 + П1 , un+1 = 62uj-1 + П2 , где vA tr_]/2 i    vA tr                                 vA tri/2

ai = —i-12-²; bi = —i+r-; c = ai + bi+1; 61 = —2—^1/2-----;

riAr²        riAr²                   0,5 r^A r + vAtr^

0,5 r/ ₂Aru 0 + 0,5 AP^j+1r12_/2Ar At1 pl ^{0,5 r}1²/2^{Ar + vAtr}1/2

П2 =

0,5Ar(Qj⁺¹- Q^j)/n + 0,5(R²- rn—_1/2)Arun - 0,5AP^j+1rn—_1/2ArAt/pl 0,5(R²- r_n²-1/2)Ar + VAtrn-1/2

62 =

___________^VAtrn-1/2___________

0,5(R²- r_n²-1/2)Ar + VAtr_n-1/2

Разностная задача (7)-( 10) при каждом фиксированном значении j, j = 1, m -1 представляет собой систему линейных алгебраических уравнений с трехдиагональной матрицей, в которой в качестве неизвестных выступают приближенные значения искомой функций u(x, t) в узлах разностной сетки, т. е. uj⁺¹, i = 0, n, j = 1, m -1.

Для решения таких систем можно использовать алгоритм Томаса, представляющий собой экономичный вариант алгоритма для метода последовательного исключения неизвестных Гаусса (метод прогонки) [18]. Согласно алгоритму Томаса решение системы (7)-(10) при каждом фиксированном значении j, j = 1, m -1 представляется в виде u+1 = a+u :1 + Д+1, i = n-1,n-2,...,1,0, un+ = 6в + П2 ,

Таблица 1

Предложенный вычислительный алгоритм также опробован на данных модели стационарно-

го течения вязкой несжимаемой жидкости в трубопроводе и d ( du)   AP _

• ——I r   I=—, 0 < ^r< R,

r dr ^ dr )     l

R

• — r=o = 0, V^Truudr = Q, dr         ₀

где Q и AP являются постоянными. Данная задача имеет точное решение

, . Q     AP Г _Р2 - 2х

u^(r) =---2+ ^⁽R - ²r ).

nR²8ц1

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

Таблица2

р = 1000 кг/м3; ц = 10 3 Паю; Q = 1,885 м3/с; AP/l = 0,2 -10 4 Па/м; R = 0,6 м, представлены в табл. 2; в ней и - вычисленные значения скорости течения при выходе на стационар по предложенному алгоритму, м/c; ust - вычисленные значения скорости течения по формуле (11), м/с.

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

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

• 1.    Лурье, М.В. Математическое моделирование процессов трубопроводного транспорта нефти, нефтепродуктов и газа / М.В. Лурье. - М.: ФГУП Изд-во «Нефть и газ» РГУ нефти и газа им. И.М. Губкина, 2003. - 336 с.

• 2.    Леонов, Е.Г. Гидроаэромеханика в бурении / Е.Г. Леонов, В.И. Исаев. - М.: Недра, 1987. -304 с.

• 3.    Басниев, К.С. Нефтегазовая гидромеханика / К.С. Басниев, Н.М. Дмитриев, Г.Д. Розенберг. - Москва-Ижевск: Институт компьютерных исследований, 2005. - 544 с.

• 4.    Лойцянский, Л.Г. Механика жидкости и газа / Л.Г. Лойцянский. - М.: Наука, 1987. - 840 с.

• 5.    Рабинович, Е.З. Гидравлика / Е.З. Рабинович. - М.: Недра, 1980. - 278 с.

• 6.    Шлихтинг, Г. Теория пограничного слоя / Г. Шлихтинг. - М: Наука, 1969. - 711 с.

• 7.    Neto, C. Boundary slip in Newtonian liquids: a review of experimental studies / C. Neto, D. Evans, E. Bonaccurso // Reports on Progress in Physics. – 2005. – Vol. 68, no. 12. – P. 2859–2897.

• 8.    Lauga, E. Microfluidics: the no-slip boundary condition in Handbook of Experimental Fluid Dynamics / E. Lauga, M.P. Brenner, H.A. Stone. – New York: Springer, 2006. – P. 1219–1240.

• 9.    Янков, В.И. Переработка волокнообразующих полимеров. Основы реологии полимеров и течение полимеров в каналах / В.И. Янков. - Москва-Ижевск: НИЦ «Регулярная и хаотическая динамика», Институт компьютерных исследований, 2008. - 264 с.

• 10.    Rao, I.J. The effect of the slip boundary condition on the flow of fluids in a channel / Rao I.J., K.R. Rajagopal // Acta Mechanica. – 1999. – Vol. 135. – P. 113–126.

• 11.    Борзенко, Е.И. Исследование явления проскальзывания в случае течения вязкой жидкости в изогнутом канале / Е.И. Борзенко, О.А. Дьякова, Г.Р. Шрагер // Вестник Томского Государственного Университета. Математика и Механика. - 2014. - № 2(28). - С. 35-44.

• 12.    Volker, J. Slip with friction and penetration with resistance boundary conditions for the Navier– Stokes equation – numerical tests and aspects of the implementation / J. Volker // J. Computational and Applied Mechanics. – 2002. – Vol. 147, Issue 2. – P. 287–300.

• 13.    Cannon, J.R. The solution of heat equation subject to the specification of energy / J.R. Cannon // Quart. Appl. Math. - 1963. - Vol. 21, № 2. - P. 155-160.

• 14.    Ионкин, Н.И. Решение одной краевой задачи теории теплопроводности с неклассическим краевым условием / Н.И. Ионкин // Дифференциальные уравнения. - 1977. - Т. 13, № 2. - С. 294304.

• 15.    Самарский, А.А. О некоторых проблемах современной теории дифференциальных уравнений / А.А. Самарский // Дифференциальные уравнения. - 1980. - Т. 16, № 11. - С. 1925-1935.

• 16.    Нахушева, В.А. Дифференциальные уравнения математических моделей нелокальных процессов / В.А. Нахушева. - М.: Наука, 2006. - 173 с.

• 17.    Самарский, А.А. Теория разностных схем / А.А. Самарский. - М. Наука, 1989. - 616 с.

• 18.    Самарский, А.А. Численные методы решения обратных задач математической физики / А.А. Самарский, П.Н. Вабищевич. - М.: Издательство ЛКИ, 2009. - 480 с.