Известия РАН. Механика жидкости и газа, 2021, № 1, стр. 3-11

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

Е. И. Борзенко a, И. А. Рыльцев a*, Г. Р. Шрагер a

a Томский государственный университет
Томск, Россия

* E-mail: ryltsev_i@ftf.tsu.ru

Поступила в редакцию 18.03.2020
После доработки 21.06.2020
Принята к публикации 21.06.2020

Полный текст (PDF)

Аннотация

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

Ключевые слова: ламинарное течение, труба с сужением /расширением, степенная жидкость, модель Оствальда-де Ваале, уравнение Пуассона для давления, преобразование координат, местное гидравлическое сопротивление

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

Для оценки влияния препятствия на распределения кинематических и динамических характеристик течения проведены многочисленные экспериментальные и численные исследования потоков в каналах и трубах с различной формой препятствия. Геометрия препятствия, рассматриваемая в данной работе, часто применяется при моделировании течений крови в стенозированных сосудах. Исследованию течений крови в сосудах с использованием ньютоновской модели реологического поведения посвящено большое количество работ [115]. Авторы экспериментальной работы [1] изучили течение жидкости в трубе для широкого диапазона чисел Рейнольдса и геометрических параметров препятствия. По результатам этих экспериментов удалось установить критические значения числа Рейнольдса для заданных геометрий, при которых происходит переход режима течения от ламинарного к турбулентному, описать структуру потока и определить потери давления, связанные с наличием препятствия в трубе. Известны исследования, например [47], где численно моделируется стационарный поток крови через стенозированную артерию. В [4] форма стеноза задается косинусоидальной функцией, а в [5, 6] – гауссовой функцией. Работа [8] посвящена изучению влияния формы стеноза на структуру течения. Авторы рассматривают четыре формы стеноза: прямоугольную, трапецеидальную и стеноз, задаваемый функцией косинуса и Гаусса. В [9] исследуются течения в каналах с несимметричными конфигурациями препятствия. В работах [10, 11] моделируются течения жидкости через два последовательно расположенных сужения, изучено влияние наличия второго сужения на распределения характеристик потока. В [1214] рассматривается нестационарное течение ньютоновской жидкости в канале с препятствием заданной формы.

Экспериментальные исследования показали, что кровь при низких скоростях сдвига проявляет неньютоновские свойства. Имеется значительно меньше работ по моделированию течений в каналах с препятствием, где в качестве исследуемой жидкости рассматривается реологически сложная среда [1621]. В [16] численно исследуется течение в крупной артерии, а в качестве рассматриваемой среды используется вязкоупругая жидкость. В [17] для описания реологических свойств течения жидкости в плоском канале авторы используют степенную модель. В [1820] исследуются пульсирующие течения неньютоновской среды, для описания которой применяются различные реологические модели, в частности, в [18] используется модель Балкли–Гершеля. В отмеченных работах по численному моделированию течений используется двумерная постановка задачи. Развитие средств численного моделирования позволило авторам [22] рассмотреть трехмерное течение в круглой трубе с препятствием, используя пакет прикладных программ ANSYS.

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

1. ПОСТАНОВКА ЗАДАЧИ

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

Рис. 1.

Область течения.

Для математического описания течения используется постановка в переменных функция тока – вихрь, которая в безразмерной форме в цилиндрической системе координат записывается в виде [23]

(1.1)
$\frac{{\partial (u{\omega })}}{{\partial z}} + \frac{{\partial ({v}{\omega })}}{{\partial r}} = \frac{B}{{\operatorname{Re} }}\left( {{{\nabla }^{2}}{\omega } - \frac{{\omega }}{{{{r}^{2}}}}} \right) + \frac{S}{{\operatorname{Re} }}$
(1.2)
${{\nabla }^{2}}{\psi } - \frac{2}{r}\frac{{\partial \psi }}{{\partial r}} = - r{\omega }$
${{\nabla }^{2}} = \frac{{{{\partial }^{2}}}}{{\partial {{z}^{2}}}} + \frac{{{{\partial }^{2}}}}{{\partial {{r}^{2}}}} + \frac{1}{r}\frac{\partial }{{\partial r}}$
$S = 2\frac{{{{\partial }^{2}}B}}{{\partial z\partial r}}\left( {\frac{{\partial {v}}}{{\partial r}} - \frac{{\partial u}}{{\partial z}}} \right) + 2\frac{{\partial B}}{{\partial z}}\frac{{\partial {\omega }}}{{\partial z}} + 2\frac{{\partial B}}{{\partial r}}\frac{{\partial {\omega }}}{{\partial r}} + \left( {\frac{{{{\partial }^{2}}B}}{{\partial {{z}^{2}}}} - \frac{{{{\partial }^{2}}B}}{{\partial {{r}^{2}}}}} \right)\left( {\frac{{\partial u}}{{\partial r}} + \frac{{\partial {v}}}{{\partial z}}} \right) + \frac{{\omega }}{r}\frac{{\partial B}}{{\partial r}}$
где безразмерные компоненты скорости и вихрь определяются формулами

$u = \frac{1}{r}\frac{{\partial {\psi }}}{{\partial r}},\quad {v} = - \frac{1}{r}\frac{{\partial {\psi }}}{{\partial z}},\quad {\omega } = \frac{{\partial {v}}}{{\partial z}} - \frac{{\partial u}}{{\partial r}}$

Реологические свойства среды описываются уравнением Освальда-де Ваале [24], в котором безразмерная эффективная вязкость вычисляется по формуле

(1.3)
$B = {{A}^{{m - 1}}}$

Здесь u, ${v}$ − аксиальная и радиальная компоненты скорости соответственно, Re = $({\rho }{{U}^{{2 - m}}}r_{0}^{m}){\text{/}}k$ – число Рейнольдса, ρ – плотность, U – среднерасходная скорость, $A = {{(2{{e}_{{ij}}}{{e}_{{ji}}})}^{{0.5}}}$, eij – компоненты тензора скоростей деформаций, k – консистенция жидкой среды, m – степень нелинейности жидкости. За масштабы обезразмеривания взяты следующие величины: длины – радиус трубы r0; скорости – среднерасходная скорость U.

В сечении Г1 профиль скорости соответствует установившемуся течению степенной жидкости в круглой трубе с постоянным расходом. На границе Г2 выполняются условия прилипания, на оси симметрии Г4 используются условия симметрии, в выходном сечении Г3 реализуются мягкие граничные условия. Таким образом, условия на границах в переменных функция тока – вихрь запишутся в виде

$\begin{gathered} {{{\text{Г}}}_{{\text{1}}}}{\text{:}}\;{\psi } = \int\limits_0^r {urdr} ,\quad {\omega } = - \frac{{\partial u}}{{\partial r}} \\ {{{\text{Г}}}_{{\text{2}}}}{\text{:}}\;\frac{{\partial {\psi }}}{{\partial n}} = 0,\quad {\omega } = \frac{1}{r}\frac{{{{\partial }^{2}}{\psi }}}{{\partial {{n}^{2}}}} \\ {{{\text{Г}}}_{{\text{3}}}}{\text{:}}\;\frac{{\partial {\psi }}}{{\partial z}} = 0,\quad \frac{{\partial {\omega }}}{{\partial z}} = 0 \\ {{{\text{Г}}}_{{\text{4}}}}{\text{:}}\;{\psi } = 0,\quad {\omega } = 0 \\ \end{gathered} $
где n – нормаль к границе Г2,$u = \left[ {(3m + 1){\text{/}}(m + 1)} \right](1 - {{r}^{{1/m + 1}}})$ при z = 0.

Граница Г2 описывается функцией

$f(z) = \left\{ \begin{gathered} 1 - \frac{\alpha }{2}\left( {1 - \cos \frac{{\pi (z - {{L}_{{\text{1}}}})}}{{{{Z}_{0}}}}} \right),\quad z \in \left[ {{{L}_{{\text{1}}}};{{L}_{{\text{2}}}}} \right] \hfill \\ 1,\quad z \in [0;{{L}_{{\text{1}}}}) \cup ({{L}_{{\text{2}}}};L] \hfill \\ \end{gathered} \right.$

Значения L1 и L2 задаются такими, чтобы исключить влияние препятствия на характер течения в окрестности Г1, Г3.

2. МЕТОД РАСЧЕТА

Решение дифференциальных уравнений для вихря и функции тока осуществляется численно с использованием метода установления, сущность которого заключается в преобразовании стационарной задачи в нестационарную. Для его реализации в уравнения (1.1) и (1.2) вводятся фиктивные производные по времени искомых функций ψ, ω и расчет по времени продолжается до обнуления введенных производных [25]. Для целей последующей разностной аппроксимации область течения с границей f(z) при преобразовании координат ${\xi } = z$, ${\eta } = r{\text{/}}f(z)$ трансформируется в прямоугольную.

Уравнения (1.1), (1.2) с введенными производными по времени в системе координат (ξ, η) имеют вид

(2.1)
$\begin{gathered} \frac{{\partial {\omega }}}{{\partial t}} + \frac{{\partial \left( {u{\omega }} \right)}}{{\partial {\xi }}} - g\frac{{\partial \left( {u{\omega }} \right)}}{{\partial {\eta }}} + \frac{1}{f}\frac{{\partial \left( {{v}{\omega }} \right)}}{{\partial {\eta }}} = \\ = \;\frac{B}{{\operatorname{Re} }}\left( {\frac{{{{\partial }^{2}}{\omega }}}{{\partial {{{\xi }}^{2}}}} + H\frac{{\partial {\omega }}}{{\partial {\eta }}} - 2g\frac{{{{\partial }^{2}}{\omega }}}{{\partial {\xi }\partial {\eta }}} + G\frac{{{{\partial }^{2}}{\omega }}}{{\partial {{{\eta }}^{2}}}} + \frac{1}{{{{f}^{2}}{\eta }}}\frac{{\partial {\omega }}}{{\partial {\eta }}} - \frac{{\omega }}{{{{f}^{2}}{{{\eta }}^{2}}}}} \right) + \frac{S}{{\operatorname{Re} }} \\ \end{gathered} $
(2.2)
$\frac{{\partial {\psi }}}{{\partial t}} + \frac{{{{\partial }^{2}}{\psi }}}{{\partial {{{\xi }}^{2}}}} + \left( {H - \frac{1}{{{{f}^{2}}{\eta }}}} \right)\frac{{\partial {\psi }}}{{\partial {\eta }}} - 2g\frac{{{{\partial }^{2}}{\psi }}}{{\partial {\xi }\partial {\eta }}} + G\frac{{{{\partial }^{2}}{\psi }}}{{\partial {{{\eta }}^{2}}}} = - f{\eta \omega }$
где $G = \frac{1}{{{{f}^{2}}}} + {{\left( {{\eta }\frac{{f{\text{'}}}}{f}} \right)}^{2}}$, $H = {\eta }\left( {2\frac{{{{f}^{{'2}}}}}{{{{f}^{2}}}} - \frac{{f{\text{''}}}}{f}} \right)$, $g = {\eta }\frac{{f{\text{'}}}}{f}$, $S = {{S}_{1}} + {{S}_{2}} + {{S}_{3}} + {{S}_{4}}$,

$\begin{gathered} {{S}_{1}} = 2\left( {\frac{1}{f}\frac{{{{\partial }^{2}}B}}{{\partial {\xi }\partial {\eta }}} - \frac{g}{{f{\eta }}}\frac{{\partial B}}{{\partial {\eta }}} - \frac{g}{f}\frac{{{{\partial }^{2}}B}}{{\partial {{{\eta }}^{2}}}}} \right)\left( {\frac{1}{f}\frac{{\partial {v}}}{{\partial {\eta }}} - \frac{{\partial u}}{{\partial {\xi }}} + g\frac{{\partial u}}{{\partial {\eta }}}} \right), \\ {{S}_{2}} = 2\left( {\frac{{\partial B}}{{\partial {\xi }}} - g\frac{{\partial B}}{{\partial {\eta }}}} \right)\left( {\frac{{\partial {\omega }}}{{\partial {\xi }}} - g\frac{{\partial {\omega }}}{{\partial {\eta }}}} \right) + 2\frac{1}{{{{f}^{2}}}}\frac{{\partial B}}{{\partial {\eta }}}\frac{{\partial {\omega }}}{{\partial {\eta }}}, \\ {{S}_{3}} = \left( {\frac{{{{\partial }^{2}}B}}{{\partial {{{\xi }}^{2}}}} + H\frac{{\partial B}}{{\partial {\eta }}} - 2g\frac{{{{\partial }^{2}}B}}{{\partial {\xi }\partial {\eta }}} + \left( {{{g}^{2}} - \frac{1}{{{{f}^{2}}}}} \right)\frac{{{{\partial }^{2}}B}}{{\partial {{{\eta }}^{2}}}}} \right)\left( {\frac{1}{f}\frac{{\partial u}}{{\partial {\eta }}} + \frac{{\partial {v}}}{{\partial {\xi }}} - g\frac{{\partial {v}}}{{\partial {\eta }}}} \right),\quad {{S}_{4}} = \frac{{\omega }}{{{{f}^{2}}{\eta }}}\frac{{\partial B}}{{\partial {\eta }}} \\ \end{gathered} $

Область решения покрывается квадратной сеткой с шагом h. Для разностного представления уравнений в преобразованной системе координат (2.1) и (2.2) применяется явная разностная схема, конвективные слагаемые в основных уравнениях аппроксимируются разностями против потока. Условия сходимости итерационного процесса записываются в виде

$\mathop {\max }\limits_{i,j} \left| {1 - \frac{{{\omega }_{{i,j}}^{{l + 1}}}}{{{\omega }_{{i,j}}^{l}}}} \right| < {\varepsilon ,}\quad \mathop {\max }\limits_{i,j} \left| {1 - \frac{{{\psi }_{{i,j}}^{{l + 1}}}}{{{\psi }_{{i,j}}^{l}}}} \right| < {\varepsilon }$
где ε – заданный итерационный параметр, l – номер шага по фиктивному времени. Значение параметра ε принималось равным 10–5. Для расчета граничных значений вихря на твердой стенке используется условие Тома [23].

Для подтверждения аппроксимационной сходимости выполнены расчеты на последовательности вложенных сеток. Рисунок 2 демонстрирует профили аксиальной скорости на последовательности сеток в сечении 2 (рис. 1). В табл. 1 приведены значения аксиальной скорости на оси симметрии в сечении Г3 для разных значений h. Результаты расчетов, представленные на рис. 2 и в табл. 1, демонстрируют аппроксимационную сходимость предложенного алгоритма. Все дальнейшие расчеты проводились на сетке с шагом h = 0.025.

Рис. 2.

Профиль аксиальной скорости Re = 1, m = 0.6, α = 0.5, Z0 = 1; 14h = 0.1,0.05, 0.025, 0.0125.

Таблица 1.

Зависимость umax от шага сетки h: Re = 1, m = 0.6, α = 0.5, Z0 = 1

h umax
0.1 1.7489
0.05 1.7518
0.025 1.7516
0.0125 1.7516
Аналитическое значение 1.75

3. ВОССТАНОВЛЕНИЕ ДАВЛЕНИЯ

Для определения поля давления используется уравнение Пуассона, которое получается в результате математических преобразований из уравнений Навье–Стокса с учетом уравнения неразрывности [23], и в физических переменных (u, ${v}$, p) в безразмерной форме принимает вид

(3.1)
$\begin{gathered} {{\nabla }^{2}}p = 2\left( {\frac{{\partial u}}{{\partial z}}\frac{{\partial {v}}}{{\partial r}} - \frac{{\partial u}}{{\partial r}}\frac{{\partial {v}}}{{\partial z}} - \frac{{{{{v}}^{2}}}}{{{{r}^{2}}}}} \right) + \\ + \;\frac{2}{{\operatorname{Re} }}\left( {\frac{{\partial B}}{{\partial z}} \times {{\nabla }^{2}}u + \frac{{\partial B}}{{\partial r}} \times {{\nabla }^{2}}{v} + \frac{{{{\partial }^{2}}B}}{{\partial z\partial r}}\left( {\frac{{\partial u}}{{\partial r}} + \frac{{\partial {v}}}{{\partial z}}} \right) + \frac{{{{\partial }^{2}}B}}{{\partial {{z}^{2}}}}\frac{{\partial u}}{{\partial z}} + \frac{{{{\partial }^{2}}B}}{{\partial {{r}^{2}}}}\frac{{\partial {v}}}{{\partial r}}} \right) \\ \end{gathered} $

В качестве масштаба давления используется величина ρU2. Уравнение (3.1) в системе координат (ξ, η) принимает вид

(3.2)
$\begin{gathered} \left( {\frac{{{{\partial }^{2}}P}}{{\partial {{{\xi }}^{2}}}} - 2g\frac{{{{\partial }^{2}}P}}{{\partial {\xi }\partial {\eta }}} + \left( {H + \frac{1}{{{{f}^{2}}{\eta }}}} \right)\frac{{\partial P}}{{\partial {\eta }}} + G\frac{{{{\partial }^{2}}P}}{{\partial {{{\eta }}^{2}}}}} \right) = 2\left( {{{S}_{5}} \times {{S}_{6}} - {{S}_{7}} \times {{S}_{8}} - {{{\left( {\frac{{v}}{{f{\eta }}}} \right)}}^{2}}} \right) + \\ + \;\frac{2}{{\operatorname{Re} }}\left( {{{S}_{9}}\left( {{{S}_{{10}}} + {{S}_{{11}}} + \frac{{{{S}_{8}}}}{{f{\eta }}}} \right) + {{S}_{{12}}}\left( {{{S}_{{13}}} + {{S}_{{14}}} + \frac{{{{S}_{6}}}}{{f{\eta }}}} \right) + {{S}_{{15}}}\left( {{{S}_{8}} + {{S}_{7}}} \right) + {{S}_{{16}}} \times {{S}_{5}} + {{S}_{{17}}} \times {{S}_{6}}} \right) \\ \end{gathered} $
где

$\begin{gathered} {{S}_{5}} = \frac{{\partial u}}{{\partial {\xi }}} - g\frac{{\partial u}}{{\partial {\eta }}}{\text{,}}\quad {{S}_{6}} = \frac{1}{f}\frac{{\partial {v}}}{{\partial {\eta }}},\quad {{S}_{7}} = \frac{{\partial {v}}}{{\partial {\xi }}} - g\frac{{\partial {v}}}{{\partial {\eta }}},\quad {{S}_{8}} = \frac{1}{f}\frac{{\partial u}}{{\partial {\eta }}} \\ {{S}_{9}} = \frac{{\partial B}}{{\partial {\xi }}} - g\frac{{\partial B}}{{\partial {\eta }}},\quad {{S}_{{10}}} = \frac{{{{\partial }^{2}}u}}{{\partial {{{\xi }}^{2}}}} - 2g\frac{{{{\partial }^{2}}u}}{{\partial {\xi }\partial {\eta }}} + H\frac{{\partial u}}{{\partial {\eta }}} + {{g}^{2}}\frac{{{{\partial }^{2}}u}}{{\partial {{{\eta }}^{2}}}} \\ {{S}_{{11}}} = \frac{1}{{{{f}^{2}}}}\frac{{{{\partial }^{2}}u}}{{\partial {{{\eta }}^{2}}}},\quad {{S}_{{12}}} = \frac{1}{f}\frac{{\partial B}}{{\partial {\eta }}},\quad {{S}_{{13}}} = \frac{{{{\partial }^{2}}{v}}}{{\partial {{{\xi }}^{2}}}} - 2g\frac{{{{\partial }^{2}}{v}}}{{\partial {\xi }\partial {\eta }}} + H\frac{{\partial {v}}}{{\partial {\eta }}} + {{g}^{2}}\frac{{{{\partial }^{2}}{v}}}{{\partial {{{\eta }}^{2}}}} \\ \end{gathered} $
$\begin{gathered} {{S}_{{14}}} = \frac{1}{{{{f}^{2}}}}\frac{{{{\partial }^{2}}{v}}}{{\partial {{{\eta }}^{2}}}},\quad {{S}_{{15}}} = \frac{1}{f}\frac{{{{\partial }^{2}}B}}{{\partial {\xi }\partial {\eta }}} - \frac{g}{{f{\eta }}}\frac{{\partial B}}{{\partial {\eta }}} - \frac{g}{f}\frac{{{{\partial }^{2}}B}}{{\partial {{{\eta }}^{2}}}} \\ {{S}_{{16}}} = \frac{{{{\partial }^{2}}B}}{{\partial {{{\xi }}^{2}}}} - 2g\frac{{{{\partial }^{2}}B}}{{\partial {\xi }\partial {\eta }}} + H\frac{{\partial B}}{{\partial {\eta }}} + {{g}^{2}}\frac{{{{\partial }^{2}}B}}{{\partial {{{\eta }}^{2}}}},\quad {{S}_{{17}}} = \frac{1}{{{{f}^{2}}}}\frac{{{{\partial }^{2}}B}}{{\partial {{{\eta }}^{2}}}} \\ \end{gathered} $

Для записи граничного условия на границе Г2 используются уравнения движения [23], которые с учетом условий прилипания на твердой стенке при η = 1 в системе координат (ξ, η) принимают вид

(3.3)
$\frac{{\partial P}}{{\partial {\xi }}} - g\frac{{\partial P}}{{\partial {\eta }}} = - \frac{1}{{\operatorname{Re} }}\left( {\frac{B}{f}\left( {\frac{{\partial {\omega }}}{{\partial {\eta }}} + \frac{{\omega }}{{\eta }}} \right) + \frac{{\omega }}{f}\frac{{\partial B}}{{\partial {\eta }}}} \right)$
(3.4)
$\frac{{\partial P}}{{\partial {\eta }}} = \frac{f}{{\operatorname{Re} }}\left( {B\left( {\frac{{\partial {\omega }}}{{\partial {\xi }}} - g\frac{{\partial {\omega }}}{{\partial {\eta }}}} \right) - {\omega }\left( {\frac{{\partial B}}{{\partial {\xi }}} - g\frac{{\partial B}}{{\partial {\eta }}}} \right)} \right)$

Из уравнений (3.3) и (3.4) получаем уравнение для расчета давления вдоль границы Г2

(3.5)
$\frac{{\partial P}}{{\partial {\xi }}} = \frac{1}{{\operatorname{Re} }}\left( {gf\left( {B\left( {\frac{{\partial {\omega }}}{{\partial {\xi }}} - g\frac{{\partial {\omega }}}{{\partial {\eta }}}} \right) - {\omega }\left( {\frac{{\partial B}}{{\partial {\xi }}} - g\frac{{\partial B}}{{\partial {\eta }}}} \right)} \right) - \left( {\frac{B}{f}\left( {\frac{{\partial {\omega }}}{{\partial {\eta }}} + \frac{{\omega }}{f}} \right) + \frac{{\omega }}{f}\frac{{\partial B}}{{\partial {\eta }}}} \right)} \right)$

Граничное условие в сечении Г1 для уравнения Пуассона при ξ = 0 определяется из уравнения (3.3), в выходном сечении давление задается равным 0, на оси симметрии используются условия симметрии $\partial P{\text{/}}\partial {\eta = 0}$.

Уравнение Пуассона для давления решается методом установления. Все пространственные производные аппроксимируются со вторым порядком точности. Для определения сходимости процесса используется неравенство

$\mathop {\max }\limits_{i,j} \left| {1 - \frac{{P_{{i,j}}^{{l + 1}}}}{{P_{{i,j}}^{l}}}} \right| < {\varepsilon }$
где ε принимается равным 10–7.

4. РАСЧЕТ МЕСТНОГО ГИДРАВЛИЧЕСКОГО СОПРОТИВЛЕНИЯ

Перепад давления между входным и выходным сечениями трубы с препятствием можно представить в виде суммы [26]

${{p}_{{\text{1}}}} - {{p}_{2}} = \Delta p = \Delta {{p}_{{{\text{тр}}}}} + \Delta {{p}_{{\text{м}}}}$
где p1, p2 – значения давления на границах Г1 и Г2 соответственно, ∆pтр – потери на трение, ∆pм – местные потери, возникающие в окрестности препятствия. Отношение потерянного на участке полного давления к динамическому напору определяет коэффициент гидравлического сопротивления [26]

$С = \frac{{\Delta p}}{{0.5{\rho }{{U}^{2}}}} = {{C}_{{{\text{тр}}}}} + {{C}_{{\text{м}}}}$

Таким образом, коэффициент местного гидравлического сопротивления рассчитывается по формуле

${{C}_{{\text{м}}}} = \frac{{\Delta {{p}_{{\text{м}}}}}}{{0.5{\rho }{{U}^{2}}}}$

Значение ∆pм вычисляется экстраполяцией давления с участков одномерного течения перед и за препятствием в сечение 1 (рис. 1). Давление на этих участках изменяется по линейному закону.

5. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Согласно результатам [21], картины стационарного течения степенной жидкости в трубе с участком сужения/расширения характеризуются зонами одномерного течения перед и за препятствием и зоной двумерного течения в окрестности препятствия с образованием циркуляционной зоны. На рис. 3 показаны типичные распределения линий тока для псевдопластичной жидкости (m < 1), демонстрирующие появление циркуляционной зоны за препятствием с увеличением числа Рейнольдса.

Рис. 3.

Распределения линий тока при m = 0.8, α = 0.5, Z0 = 1: (a), (б) – Re = 1, 25.

На рис. 4а представлены зависимости перепада давления от числа Рейнольдса для ньютоновской жидкости. Такая область решения использовалась в работах [10, 27], где представлены результаты численного моделирования. Экспериментальное исследование течения в трубе данной геометрии описано в [1], где установлено, что значение числа Рейнольдса выше 250 при заданных геометрических параметрах приводит к потере устойчивости ламинарного режима. Наблюдается хорошее согласование результатов расчетов с данными, представленными в [1, 10, 27].

Рис. 4.

(а) Зависимость перепада давления от Re при m = 1, α = 0.333, Z0 = 4, L = 32, L1 = 12: 1 – настоящая работа; 24 – исследования [1, 10, 27]; (б) распределение давления вдоль оси z при Re = 12.5, α = 0.5, Z0 = 2: 58m = 0.8, 1, 1 ([11]), 1.2.

На рис. 4б показаны изменения давления вдоль оси симметрии трубы для трех значений показателя нелинейности. Полученное распределение давления при m = 1 практически совпадает с данными из [11]. Максимальное значение перепада давления реализуется в случае дилатантной жидкости (m = 1.2). Анализ графиков, представленных на рис. 4б, позволяет сделать вывод, что для псевдопластичной жидкости (m = 0.8) перепад давления между входным и выходным сечениями трубы уменьшается по сравнению с течением ньютоновской жидкости. Такое поведение распределения давления согласуется с изменениями распределений эффективной вязкости в зависимости от параметра нелинейности, представленными в [21].

На рис. 5 показаны распределения давления для псевдопластичной жидкости при трех значениях параметра α. Перепад давления в областях одномерного течения постоянный и соответствует течению Пуазейля в круглой трубе. Увеличение глубины перекрытия способствует росту зоны двумерного течения за препятствием. Формирование циркуляционной зоны за препятствием [21] приводит к возникновению области с минимальным давлением, локализованной в окрестности максимального сужения трубы.

Рис. 5.

Поле давления Re = 10, m = 0.8, Z0 = 1: (а–в) – α = 0.15, 0.333, 0.5.

Влияние параметра нелинейности для двух значений Re на коэффициент местного гидравлического сопротивления представлено на рис. 6. С уменьшением m коэффициент местного гидравлического сопротивления падает, что обусловлено ростом среднего уровня вязкости в области двумерного течения, влияние показателя нелинейности усиливается с уменьшением числа Рейнольдса.

Рис. 6.

Зависимости коэффициента местного гидравлического сопротивления от показателя нелинейности α = 0.5, Z0 = 1: 1, 2 – Re = 1, 10.

Зависимости См от числа Рейнольдса для разных значений m и геометрических параметров α, Z0 представлены на рис. 7. Графические зависимости показывают, что с ростом числа Рейнольдса коэффициент местного гидравлического сопротивления монотонно уменьшается. Уменьшение глубины перекрытия канала α также приводит к уменьшению значений коэффициента местного гидравлического сопротивления. Влияние параметра Z0 становится существенным при малых значениях α.

Рис. 7.

Зависимость См от Re: (а–в) – m = 0.8, 1, 1.2; 16 – (α, Z0) = (0.5, 2), (0.5, 1), (0.333, 2), (0.333, 1), (0.15, 2), (0.15, 1).

ЗАКЛЮЧЕНИЕ

Выполнено математическое моделирование течения степенной жидкости в трубе с препятствием. Разработан и реализован численный алгоритм решения поставленной задачи. Проведены параметрические исследования динамических характеристик потока для показателя нелинейности $0.6 \leqslant m \leqslant 1.4$ и числа Рейнольдса в диапазоне от 1 до 100. Вычислены значения коэффициента местного гидравлического сопротивления и выявлены тенденции его изменения при исследовании основных параметров задачи. Установлено, что с увеличением числа Рейнольдса коэффициент сопротивления уменьшается. Рост показателя нелинейности способствует увеличению коэффициента cопротивления. Уменьшение глубины перекрытия канала приводит к уменьшению значения коэффициента.

Исследование выполнено за счет гранта РНФ (проект № 18-19-00021).

Список литературы

  1. Young D.F., Tsai F.Y. Flow characteristics in model of arterial stenoses-I. Steady flow // J. Biomech. 1973. V. 6. № 4. P. 395–402.

  2. Forrester J.H., Yang D.F. Flow through a converging-diverging tube and its implications in occlusive vascular disease-I // J. Biomech. 1970. V. 3. № 3. P. 297–305.

  3. Liepsch D. , Singh M., Lee M. Experimental analysis of the influence of stenotic geometry on steady flow // Biorheology. 1992. V. 29. № 4. P. 419–431. https://www.scopus.com/authid/detail.uri?authorId=19535741000&amp;eid=2-s2.0-0027093633

  4. Deshpande M.D., Giddens D.P., Mabon R.F. Steady laminar flow through modelled vascular stenosis // J. Biomech. 1976. V. 9. № 4. P. 165–174.

  5. Lee J.S., Fung Y.C. Flow in locally constricted tubes at low Reynolds numbers // J. Appl. Mech. 1970. V. 37. № 1. P. 9–16.

  6. Lee T.S. Steady laminar fluid flow through variable constrictions in vascular tubes // J. Fluid Eng-T ASME. 1994. V. 116. № 1. P. 66–71.

  7. Pontrelli G. Blood flow through an axisymmetric stenosis // P. I. Mech. Eng. H. 2001. V. 215. № 1. P. 1–10.

  8. Banerjee M.K., Nag D., Ganguly R., Datta A. Hemodynamics in stenosed arteries effects of stenosis shapes // Int. J. Comp. Meth. 2010. V. 7. № 3. P. 397–419.

  9. Mandal D.K., Manna N.K., Chakrabarti S. Influence of different bell-shaped stenosis on the progression of the disease atherosclerosis // J. Mech. Sci. Technol. 2011. V. 25. № 8. P. 1933–1947.

  10. Lee T.S., Liao W., Low H.T. Numerical simulation of turbulent flow through series stenosis // Int. J. Numer. Meth. Fl. 2003.V. 42. № 7. P. 717–740.

  11. Huang H., Lee T.S., Shu C. Lattice-BGK simulation of steady flow through vascular tubes with double constrictions // Int. J. Numer. Method. H. 2006. V. 16. № 2. P. 185–203.

  12. O'Brien V., Ehrlich L.W. I. Simple pulsatile flow in an artery with a constriction // J. Biomech. 1985. V. 18. № 2. P. 117–127.

  13. Tu C., Deville M., Dheur L., Vanderschuren L. Finite element simulation of pulsatile flow through arterial stenosis // J. Biomech. 1992. V. 25. № 10. P. 1141–1152.

  14. Paul M.C., Molla M.M. Investigation of physiological pulsatile flow in a model arterial stenosis using large-eddy and direct numerical simulations // Appl. Math. Mod. 2012. V. 36. № 9. P. 4393–4413.

  15. Егоров В.А., Регирер С.А., Шадрина Н.Х. Особенности пульсирующего течения крови через резистивные кровеносные сосуды // Изв. РАН. МЖГ. 1994. № 2. С. 83.

  16. Leuprecht A., Perktold K. Computer Simulation of Non-Newtonian Effects on Blood Flow in Large Arteries // Comput. Method Biomec. 2001. V. 4. № 2. P. 149–163.

  17. Gao W., Liu R., Duan Y. Numerical investigation on non-newtonian flows through double constrictions by an unstructured finite volume method // J. Hydrodyn. 2009. V. 21. № 5. P. 622–632.

  18. Tu C., Deville M. Pulsatile flow of non-Newtonian fluids through arterial stenosis // J. Biomech. 1996. V. 29. № 7. P. 899–908.

  19. Mukhopadhyay S., Mandal M.S., Mukhopadhyay S. Effects of variable viscosity on pulsatile flow of blood in a tapered stenotic flexible artery // Math. Method Appl. Sci. 2019. V. 42. № 2. P. 488–504.

  20. Achab L. Numerical simulations of the pulsatile blood flow in narrowing small vessels using different rheological models // J. Phys. Conf. Ser. 2019. V. 1294. № 2.

  21. Рыльцев И.А., Рыльцева К.Е., Шрагер Г.Р. Кинематика течения степенной жидкости в трубе переменного сечения // Вестн. Том. гос. ун-та. Математика и механика. 2020. № 63. С. 125–138.

  22. Jung H., Choi J.W., Park C.G. Asymmetric flows of non-Newtonian fluids in symmetric stenosed artery // Korea-Aust. Rheol. J. 2004. V. 16. № 2. P. 101–108.

  23. Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.

  24. Ostwald W. Ueber die rechnerische Darstellung des Strukturgebietes der Viskosität // Kolloid Zeitschrift. 1929. V. 47. № 2. P. 176–187.

  25. Годунов С.К., Рябенький В.С. Введение в теорию разностных схем. М.: Физматгиз, 1962. 340 с.

  26. Идельчик И.Е. Справочник по гидравлическим сопротивлениям / под ред. М.О. Штейнберга. 3-е изд. М.: Машиностроение, 1992. 672 с.

  27. Zendehbudi G.R., Moayeri M.S. Comparison of physiological and simple pulsatile flows through stenosed arteries // J. Biomech. 1999. V. 32. № 9. P. 959–695.

Дополнительные материалы отсутствуют.