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

УДК 531.7

ТЕЧЕНИЕ НЕНЬЮТОНОВСКОЙ ЖИДКОСТИ В КАНАЛЕ ПРЯМОУГОЛЬНОГО СЕЧЕНИЯ

В. Ф. Никитин

Федеральный научный центр Научно-исследовательский институт системных исследований

Российской академии наук, Московский государственный университет им. М. В. Ломоносова, i_nikitin@list.ru

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

Ключевые слова: дилатантная жидкость, канальная модель, пористая среда.

NON-NEWTONIAN FLOW IN RECTANGULAR CROSS-SECTION CHANNEL

V. F. Nikitin

System Research Institute, Russian Academy of Sciences, Moscow State University, i_nikitin@list.ru

The paper investigates stationary dilatant fluid flow with a power function for viscosity in a rectangular cross section pipe driven by a specified pressure difference. The numeric study has generated an interpolation equation for the flow rate vs. other parameters.

Keywords: dilatant fluid, channel model, porous medium.

Введение

Проблема упрощения описания гидродинамики и прочих свойств пористой среды существовала с XIX века, когда была предложена модель Дарси, не детализирующая среду. В первой трети XX века была предложена модель Кармана — Козени [1-2], описывающая пористую среду как засыпку модельных шаров. Дальнейшие исследования показали, что модель Кармана — Козени не всегда адекватно описывает течение в пористой среде, и был предложен подход к описанию пористой среды как системы каналов [3-4].

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

В дальнейшем канальный подход оказался достаточно перспективен, но сложность устройства порового пространства приводит к необходимости постоянного совершенствования моделей подобного класса. Кроме того, встает проблема связи канальной модели с реальной картиной порового пространства, возможное решение которой описано, но оно годится далеко не для всех сред, поскольку форма поперечного сечения в большинстве таких канальных моделей фиксирована [6-7]. Вместе с тем важную проблему представляет собой моделирование резервуаров трудноизвлекаемой нефти, в которых зачастую поры представляют собой систему плоских микротрещин [8-10]. В таких случаях нельзя рассматривать поперечное сечение канала круглым, квадратным или звездчатым [5].

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

плане, чем предварительный расчет для сечений различной формы и жидкости с различными реологическими свойствами с дальнейшим получением алгебраической зависимости между потоком в канале и продольным градиентом давления. Подобная задача, впрочем, тоже может содержать большое число параметров, и поэтому следует их разумно ограничить. Например, форму поперечного сечения канала можно выбрать прямоугольной (что описывает как простые каналы квадратного сечения, так и узкие трещины). Жидкость же возьмем степенную, зависимость вязких напряжений в которой от скорости деформации выражается степенным законом, известным иначе как закон Оствальда — де Вела [8-9]. Этот закон единообразно описывает ньютоновскую, дилатантную и псевдопластическую жидкости. Это позволяет ограничить число параметров реологии.

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

Реологическая модель

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

тИ = Це^Ц)^, (1)

- его нор-ормации, то

где Т11 — тензор вязких напряжений, — девиатор тензора скоростей деформации ма, // — коэффициент динамической вязкости. Если вязкость не зависит от скорости дес жидкость называется линейно-вязкой, или ньютоновской.

Для несжимаемой жидкости девиатор тензора вязких напряжений совпадает с самим тензором, и его первый инвариант равен нулю. Норма тензора строится на основе второго инварианта, который имеет вид е'^е'ц. Будем использовать следующий вид нормы [9]:

= 7. (2)

Степенной закон для вязкости обобщает модели псевдопластической, ньютоновской и дилатантной жидкостей:

= , (3)

при этом степень а определяет вид жидкости, 0 < а < 1 — псевдопластическая, а =1 — ньютоновская, а > 1 — дилатантная. Первая характеризуется падением вязкости с возрастанием скорости деформации (например, сдвига) — пример среди пищевых продуктов кетчуп. Вторая имеет постоянную вязкость (обычная жидкость). Для дилатантной жидкости характерно возрастание вязкости с увеличением скорости деформации, примером является суспензия кукурузной сечки или болотная трясина. В формуле (3) /.¿о имеет размерность обычной вязкости и равно вязкости жидкости при 7 = 70. Разделение этих

параметров сделано для облегчения обезразмеривания, поскольку реально вязкость зависит не от трех

1 —а

независимых параметров, а от двух — степени и комплекса 0 0 . Постановка задачи

Поперечное сечение канала прямоугольное, расположено в плоскости 0у2, причем

0< у< 2г, 0< г< 2Хг, А<1. (4)

При А = 1 сечение квадратное, г является полушириной сечения (аналог большой полуоси). Обозначим область поперечного сечения , а ее границу .

Скорость жидкости имеет лишь одну компоненту их = и(у,г). В этом случае уравнение неразрывности выполняется автоматически, а уравнение импульса имеет следующий вид:

д ( ди\ д ( ди\ др ^ _

дУ^дУ) = аХ = -г'

и = 0, (У,г)€ дп. (5)

На стенках выполняется при этом условие прилипания. Выбор знака градиента давления сделан с тем, чтобы при положительном значении Г скорость также была положительна. Второй инвариант тензора скоростей деформации (2) имеет следующий вид:

• (6)

Для того, чтобы получить поток через поперечное сечение, следует проинтегрировать скорость и:

q = udtt. (7)

J п

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

r — половины ширины сечения («большой полуоси»);

Л — отношения длин малой стороны к большой (толщины к ширине); — заданного градиента давления в поперечном направлении;

Но — характерной вязкости;

70 — характерной скорости деформации (обратной величине к характерному времени релаксации неньютоновской жидкости, для ньютоновской жидкости этот параметр не имеет смысла);

а — показателя степени, для ньютоновской жидкости — а = 1.

Обезразмеривание задачи

Для уменьшения числа параметров, от которых зависит решение, а в идеале — для получения аналитической формулы для потока с точностью до константы, следует обезразмерить задачу. В качестве характерной шкалы расстояния будем брать r, а шкалу скорости пока обозначим U, ее значение в дальнейшем уточним. Будем обозначать безразмерные значения величин так е, как размерные, но с чертой сверху. При подстановке шкал параметров в уравнение (5) с учетом (6) и (3), будем иметь следующее уравнение:

д f_dU\df_dU\. ._ _. =

ддЫ) + (y'z)efi;

U = 0, (y, z)edñ; (8)

' du\ 2 í du\

м) V dzj

V =

2

. (9)

Для того чтобы выполнялось (8), (9), следует избрать шкалой вязкости /.¿о7о ° (П/г, а шкалу скорости связать с остальными параметрами следующим образом:

П па-1 /ГГ°

^0-;-г = Г, или и = -— • (10)

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

0< у< 2, 0< г< 2А, (11)

и тем самым решение безразмерной задачи (8), (9) оказывается зависящим лишь от двух параметров и.

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

/■2 />2А

Ц = Ц (А,а) = I I ийгйу, (12)

размерный же выразится через него как

- т,{\ '0 \ Л

00

q = Ur2q (А,«) = q (А,«) • -^ r2 =

=q ^ -)• = q № -)• г1" (13)

112 В. Ф. Никитин

Течение неньютоновской жидкости в канале прямоугольного сечения

В формуле (13) использованы различные представления потока от параметров неньютоновской жидкости и размера ширины канала. В последнем представлении специально выделена зависимость от градиента давления, которая имеет вид Г1/о, для дилатантной жидкости показатель степени при градиенте меньше 1. Можно также заметить, что степень зависимости от радиуса канала падает с 4 до 3 (в пределе) с ростом дилатансии. В предыдущем представлении выделен безразмерный множитель, позволяющий оценить отклонение для неньютоновской жидкости: первый множитель (безразмерный) определяется из решения поставленной выше безразмерной задачи, второй (с размерностью объемного потока) характерен как для неньютоновской, так и для ньютоновской жидкости, третий (безразмерный) отличен от единицы лишь для неньютоновской жидкости.

Расчетная схема и расчетные параметры

Расчетная схема задачи имеет следующие черты:

- Схема равномерна по обеим координатам, и содержит ЫуЫх ячеек, шаг сетки в каждом пространственном направлении Ну = 2/Ыу, Нх = 2 А/Мх.

- Скорость вычисляется в центрах ячеек, вязкость на гранях.

- Граничное условие Дирихле реализуется с помощью экстраполяции второго порядка, что приводит к нелинейной схеме второго порядка.

- Линеаризация нелинейной схемы проводится следующим методом: по полученному с предыдущей итерации решению рассчитывается вязкость (9), с заданным значением вязкости рассчитывается скорость (8), а затем ее интеграл (12). Это одна внешняя итерация.

- Внутренние итерации служат для расчета скорости на основе СЛАУ, полученного дискретизацией уравнения (8). Нами используется метод Б1-Св81аЬ [10], поскольку при переменной вязкости матрица системы несимметрична.

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

- Начальная скорость для самой первой итерации берется равной нулю. Начальное значение вязкости для получения решения на первой внешней итерации — единица.

- На внутренних итерациях решение ищется с точностью ец, на внешних с точностью ех. Критерий остановки внешних итераций — малая разница между начальным и конечным значением скорости во всех ячейках сетки.

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

Этим методом можно решать подобную систему при 0 < а < 2. При а и 2 ускорить простую итерацию можно методом установления. При а > 2 метод не работает (расходится на внешних итерациях), в том числе метод установления.

Аналитическое решение для малого параметра соотношения сторон

В том случае, когда параметр А мал, прямоугольное сечение представляет собой узкую щель, вытянутую вдоль координаты у. В этом случае влияние поперечных стенок этой щели на течение представляется малым по сравнению с влиянием продольных стенок, а градиенты распределенных параметров вдоль г значительно большими по сравнению с градиентами вдоль у в большинстве точек сечения. Для получения решения в пределе 0 (более точно, при 2 1) в уравнениях (8) и (9) оставляем лишь частные производные по г.

д ( дй\

д )=_1, 0 < г < 2А; и = 0, г = 0, г = 2А; (14)

я~ \ дг )

г

V =

2 ди о—

\дг) дг

(15)

Система (14)-(15) может быть решена уже аналитически, для любого значения показателя степени. Решение выглядит как:

й = - г - А 11+1), д = г - А . (16)

Интегрируя решение (16) по сечению согласно формуле (12), получим следующее выражение для

безразмерного потока д(А,а) при малом значении А:

Ц = 2^+2 (17)

Тем самым, для жидкости Ньютона безразмерный поток будет равен (4/3) А3, при а = 2 он равен (8/5) А25, а в пределе а ос он равен 2А2. Поскольку А мало, подобная зависимость говорит о возрастании этой величины при возрастании с выходом на константу, зависящую от .

Результаты расчетов: распределенные параметры

Были проведены расчеты функции Ц(А,а) для А = {0.1,0.2, ..., 0.9, 1.0}, и при каждом из этих значений А для а € {1.00, 1.02, ..., 1.98, 2.00}, функция эта имеет достаточно простую форму и к ней можно подобрать соответствующую интерполяцию. Выбор диапазона параметра обусловлен тем, что нам требуется рассчитать жидкость с дилатансией.

На рис. 1 изображены зависимости и(у ,х) (слева) и Д(у,х) (справа) при А = 1 (для квадратного сечения) для а =1.0 (жидкость Ньютона) и а = 1.5, а = 2.0 (дилатантная жидкость). Для а =1.0 вязкость не изображена, так как равна константе (единице).

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

На рис. 2 изображены те же зависимости, но для прямоугольного сечения при А = 0.5. При = 1 .0 вязкость не изображена, она равна константе (единице).

Рис. 2 показывает те же тенденции, что и рис. 1, но помимо этого видно, что область пониженного значения безразмерной вязкости не совпадает с диагоналями прямоугольника, а отходит от каждого угла под 45° внутрь сечения, в основной же части прямоугольного сечения зона пониженной вязкости занимает линию посередине сечения между пересечениями этих угловых зон. Тем самым видно, что при уменьшении параметра А следует ожидать, что в основной зоне течения параметры зависят прежде всего лишь от той координаты, по которой идет сокращение размера, влияние же узких торцов прямоугольного сечения распространяется на область, приблизительно равную размерам этих торцов.

Результаты расчетов: аналитическое решение для узкого сечения

На рис. 3 изображены графики зависимости и(г) и Д(г) для различных значений а, полученные для А = 0.1 по формулам (16). При этом значение степени дилатантной жидкости может и превышать 2, и из графиков видно, что с ростом безразмерная скорость увеличивается, а вязкость уменьшается. На середине отрезка вязкость падает до нуля, если только жидкость не ньютоновская, причем для небольшого отклонения от ньютоновской жидкости это падение к нулю наиболее резкое. Если же значение велико, то на середине отрезка профиль скорости испытывает сильный излом.

Результаты расчетов: безразмерный поток

На рис. 4 изображены зависимости Ц(а) для различных значений А в диапазоне [0.5,1.0]. Само значение показателя степени в диапазоне [1.0,2.0].

Из рис. 4 видно, что при выбранном нами методе обезразмеривания функция Ц(а) достаточно гладкая и может быть интерполирована с использованием различных методик (к примеру, полиномов). Интерполяционные коэффициенты этой функции зависят от , их также можно интерполировать, в итоге получив достаточно точную аналитическую зависимость, что и является целью этой работы. С увеличением при заданном безразмерный поток возрастает, как и с увеличением при заданном .

На рис. 5 изображены зависимости Ц(а) для различных значений А в диапазоне [0.1,0.5]. Для удобства ось ординат имеет логарифмический масштаб. Само значение показателя степени в диапазоне [ 1.0,2.0] соответствует рассчитанным численно значениям. Добавлены также пунктирные кривые, соответствующие аналитическому расчету при малых значениях .

Из рис. 5 видно, что для малых значений аналитическое решение действительно близко к численному. При этом численное решение всегда ниже аналитического. Во-вторых, зависимости численного и аналитического решения от параметра а выглядят подобными. Это может сыграть большую роль в подборе интерполяционной формулы для зависимости Ц(А,а).

1

0.9 0 8 0.7 0.6 0.5 0.4 0.3 0.2 0 1 0

1

0.9 0 8 0.7 0.6 0.5 0.4 0.3 0.2 0 1 О

и : а = 1.0; А = 1.0

0 9 0.8 0.7 0 6 0.5 0 4 0.3 0.2 0.1 О

04 0.6 0.8 1

а = 1.5; А = 1.0

0 9 0.8 0.7 0 6 0.5 0 4 0.3 0.2 0.1 О

и : а = 2.0; А = 1.0

1

0.9 О 8 0.7 0.6 0.5 0.4 О 3 0.2 0.1 О

О 9 0.8 0.7 О 6 0.5 04 0.3 0.2 0.1

0.2 0.4 0.6 0.8 1

Д : а = 1.5; А = 1.0

М|[ '■У

1

О 9 0.8 0.7 О 6 0.5 04 0.3 0.2 О 1

0.2 0.4 0.6 0.8 1

Д : а = 2.0; А = 1.0

Рис. 1. Контуры безразмерных функций скорости и вязкости в поперечном сечении потока для различных значений параметра а (ньютоновская и дилатантная жидкости). Квадратное поперечное сечение = 1 . 0

Интерполяционная формула для безразмерного потока

В основу построения формулы для безразмерного потока положим формулу (17), соответствующую решению при малых значениях А, эта малость подтверждается рис. 5. Для оценки отклонения численного решения от аналитического воспользуемся тем, что аналитическое решение при выглядит при малых А как двускатная крыша с прямыми наклонными сторонами (рис. 3, слева). Из результатов на рис. 2 видно, что область пониженного значения безразмерной вязкости не совпадает

-0.2

-0.2

-0.2

0 0.2 0.4 0.6

и : а = 1.0; А = 0.5

О 0.2 0.4 0.6

и : а = 1.5; А = 0.5

О 0.2 0.4 0.6

и : а = 2.0; А = 0.5

-0.2

-0.2

О 0.2 0.4 0.6

Д : а = 1.5; А = 0.5

О 0.2 0.4 0.6

Д : а = 2.0; А = 0.5

Рис. 2. Контуры безразмерных функций скорости и вязкости в поперечном сечении потока для различных значений параметра а. Квадратное поперечное сечение А =1.0

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

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

0.05 0.1 0.15

I

и(г) для А = 0.1

0.8

0.6

~а= 10 а = 1 5 а = 2.0 -а = 25 -а = 5 0 -а-50

/(г) для А = 0.1

Рис. 3. Графики зависимости и(г) и //(г) для различных значений а, полученные для А = 0.1 по аналитическим формулам

Рис. 4. Зависимость безразмерного потока ц(А ,а) как графики зависимости от а для различных А при 0. 5

ее вычисляется как

Ц = 1 • 2АН (2 - 2А) + 1 • 4А2Н = 2АН ^ 1-3) . (18)

Само значение высоты Н получается из аналитического решения как максимум функции и (16) при г = А. Подставив максимум при а оо, имеем Н = А. Тем самым из (18) получим следующую оценку (скорее всего, точную) потока при и произвольного значения :

? = 2А2(1-3) . (19)

Теперь воспользуемся тем, что согласно рис. 5 зависимости аналитического и численного решения от а близки к подобным, и вместо 2А2 (решения для а сю) введем в формулу (19) аналитическое реше-

Рис. 5. Зависимость безразмерного потока -(А,а) как графики зависимости от а для различных А при А < 0.5. Пунктирными кривыми соответствующего цвета показан результат аналитического решения для малого значения А

ние (17) для произвольного а. Кроме того, добавим предмет дальнейшей интерполяции — функцию, зависящую от , входящую в интерполяционное решение в составе малой поправки:

4а

-Л2+« I I — Л(п\

-ш (А,а,А(а)) =

2а + 1

(

Ж«)*) .

(20)

Для согласования (20) с решением (19) при а оо нужно, чтобы А(оо) = 1, это условие может подсказать форму интерполяционной зависимости уже для этой функции.

Расчет численного решения был проведен в узлах, которые можно расположить как элементы матрицы qi,j, i = 1,...,п0, / = 1,...,Пд; при этом соответствующие значения А,а соответствуют векторам щ, X/. Для заданного значения индекса i и соответствующего щ оптимальное значение Ai достигается минимизацией целевой функции (метод наименьших квадратов с весом) [11]:

ПА

'фi = ^2 № (А /) [-,/ - qint{Xj ,ai,A)

/=1

где № (А) — положительный вес данных, служащий для более оптимального учета разности расчетных данных и интерполяции в случае сильного влияния переменной на функцию. Поскольку зависимость qint(А) линейная, то оптимальное значение А = Ai получается по следующей формуле:

^ - 2 ± № (А/) / ^(Х/ ,«г,А)] • = 0,

дА

/=1

(л / )л /

/=1

п

3+У

-Ц

4а i

2а1 + 1 ; п

2+ — А;

2ш + 1 ;

О" а/]

= 0,

(4/3)«,- ^ 6+^ ^ з+—

А - и3^^ № (\/ )А; "i = £ Г (\/ )А; "i

2а i + 1

/=1

/=1

4а i

2ш + 1 ;

2+ —

qi,/

(21)

2

Оценить точность приближения можно как:

е = тах

Яш {А/,«г,Аг)

(22)

На рис. 6 изображена функция Л(а), рассчитанная по формуле (21) для трех различных функций веса № (А). Там же приведена точность для каждого из весов, рассчитанная по формуле (22).

Рис. 6. Зависимость А(а) в интерполяционной формуле (20), рассчитанная для интервала 1 < а < 2 для различного значения веса. В легенде приведены соответствующие значения максимального относительного отклонения интерполяционной кривой от значений, полученных численным расчетом

1

Из рис. 6 видно, что кривые А( ) для различного веса весьма близки друг к другу, но расчет максимального отклонения е по формуле (22) показывает, что лучшим приближением из протестированных является кривая, полученная с весом № = А-2, относительное отклонение между расчетным и интерполированным Я при этом составляет 3,2 %.

Далее эту кривую А (а), полученную расчетом (21), следует интерполировать, чтобы получить окончательную формулу не в табличном виде, а в виде выражения. Для этого воспользуемся тем, что теоретически А( ) = 1 , и представим интерполяцию в виде рациональной функции заданной степени 5 > 1:

1 5 —

Ам Н = 1 + -в£-|. (23)

к=1

Будем искать интерполяцию линейным методом наименьших квадратов при заданном значении модифицирующей степени В, что приведет к следующей системе линейных уравнений относительно коэффициентов —:

2 К °В 4 )+ё *) *

5 «и ] «и В I Л П

Е-Е^+т = Е^^• ^ (24)

к=1 1=1 I 1=1 1

Матрица этой системы симметрична, ее можно решить с помощью разложения Холецкого.

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

/3 = 0.014, -1= 0.79356, -2 = -0.01974,

так что окончательная интерполяционная формула будет выглядеть как

_0.014 /0.7935 0.01974 \

- /л Ч 4 а ,2+] Ф„, (Л H^+fA .

(25)

При этом максимальное относительное отклонение составляет 3,094 %, т. е. даже меньше, чем для рассчитанных с весом Лоптимальных значений А«(а). Это говорит о том, что дополнять полином (23) членами следующего порядка смысла нет.

На рис. 7 изображено значение интерполированного потока Ч« в виде карты в пространстве переменных , (слева) и относительное отклонение его от рассчитанного в процентах = \Ч — Яш\ /Я ■ 100% справа.

Interpolated flux

Relative error, %

■ 1 2

0.9 1.9

0.8 1.8

0.7 1.7

0.6 1.6

0.5 е 1.5

0.4 1.4

0.3 1.3

0.2 1.2

■ 0.1 11

qint (A,a)

0.2 0.4 0.6 0.8 \

e = \q- -qint\/q- 100%

Рис. 7. Безразмерный поток в трубе прямоугольного поперечного сечения в зависимости от соотношения сторон и показателя степени дилатантной жидкости (слева). Относительное отклонение рассчитанного и интерполированного значения потока, в процентах (справа)

Из рис. 7 видно, что относительная ошибка интерполяции максимальна для ньютоновской жидкости и соотношения сторон 0.6-0.7, а также около 1.0, но и там она не превышает 3, 1 %. В прочих областях, в том числе и при увеличении значения степени , отклонение меньше. Это говорит о достаточно хорошем качестве интерполяционной формулы (25).

Выводы

Нами получено решение нелинейной в общем случае задачи о течении несжимаемой неньютоновской степенной жидкости по каналу прямоугольного поперечного сечения под действием заданного градиента давления. Решение для величины объемного потока жидкости через поперечное сечение зависит от 2 безразмерных параметров, отношения меньшей стороны сечения к большей и показателя степени в реологическом законе а. Численное решение получено в диапазоне параметров А € [0 .1,1. 0], а € [1. 0,2 .0] (ньютоновская либо дилатантная жидкости) с интервалами 0 .1 и 0 .02 соответственно, аналитическое получено в приближении А2 << 1 для произвольного значения степени а. По рассчитанным значениям на основе аналитического решения выведена интерполяционная формула, связывающая безразмерный поток с параметрами , , содержащая только 3 эмпирических коэффициента и отклоняющаяся не более чем на 3.1% от полученных численным расчетом значений в указанном диапазоне параметров. Это решение в дальнейшем может быть использовано как одна из частей имеющей большую перспективу канальной модели пористой среды.

Работа выполнена при поддержке РФФИ (грант № 16-29-15080).

ЛИТЕРАТУРА

1. Kozeny J. Royal Academy of Science, Vienna // Proc Class I. 1927. Vol. 136. Р. 127.

2. Carmen P. C. Fluid flow through granular beds // Trans Inst Chem. Eng., London, 1937. Vol. 15. P. 150.

3. Owen J. E. The resistivity of a fluid-filled porous body // Trans. AIME. 1952. Vol. 195. P. 169.

4. Fatt I. The network model of porous media // Trans. AIME. 1956. Vol. 207. P. 144.

5. Yale D. P. Network modeling of flow, storage and deformation in porous rocks. A dissertation submitted to the department of geophysics and the committee of graduate studies of Stanford university in the partial requirements for the degree of Doctor of Philosophy. August, 1984.

6. Bakke S., 0ren P. E. 3-D pore modelling of sandstones and flow simulations in the pore networks // SPE Journal. 1997. June. V. 2.

7. Silin D., Patzek T. Pore space morphology analysis using maximal inscribed spheres // Physica. A 371. 2006. Р. 336 - 360. doi: 10.1016/j.physa. 2006.04.048.

8. Бетелин В. Б., Никитин В. Ф., Смирнов Н. Н., Михальченко Е. В., Скрылева Е. И., Стамов Л. И., Тюренкова В. В. Компьютерный керносимулятор — подходы и методы // Вестн. кибернетики. 2015. № 4. С. 33-44.

9. Никитин В. Ф., Стамов Л. И., Михальченко Е. В. Трехмерное математическое моделирование течения вязких жидкостей в многосвязной системе каналов и пор // Вестн. кибернетики. 2016. № 2. Р. 128-138.

10. Козлов И. В., Скрылева Е. И. Математическое моделирование и обработка эксперимента по вытеснению нефти водой из неокомских песчаников // Вестн. кибернетики. 2016. № 2. С. 139-146.

11. Scott Blair G. W., Hening J. C., Wagstaff A. The Flow of Cream through. Narrow Glass Tubes // J Phys Chem. 1939. № 7 (43). Р. 853-864.

12. Шульман 3. П. Беседы о реофизике. Минск, 1976. 96 с.

13. Vorst H. A., van der. Iterative Krylov Methods for Large Linear System. Cambridge University Press, 2003. 221 р.

14. Линник Ю. В. Метод наименьших квадратов и основы математико-статистической теории обработки наблюдений. 2-е изд. М. : Гос. изд. физ.-мат. лит., 1962. 354 с.