Известия РАН. Механика жидкости и газа, 2021, № 1, стр. 3-11
ПОТЕРИ ДАВЛЕНИЯ ДЛЯ ТЕЧЕНИЯ СТЕПЕННОЙ ЖИДКОСТИ В ТРУБЕ ПЕРЕМЕННОГО СЕЧЕНИЯ
Е. И. Борзенко a, И. А. Рыльцев a, *, Г. Р. Шрагер a
a Томский государственный университет
Томск, Россия
* E-mail: ryltsev_i@ftf.tsu.ru
Поступила в редакцию 18.03.2020
После доработки 21.06.2020
Принята к публикации 21.06.2020
Аннотация
Выполнено численное моделирование установившегося ламинарного течения несжимаемой степенной жидкости в трубе с препятствием заданной формы. Для математического описания процесса используются уравнение переноса вихря и уравнение Пуассона для функции тока, при этом реологические свойства среды описываются степенным законом Оствальда-де Ваале. Получение стационарного решения сформулированной задачи осуществляется методом установления на основе конечно-разностной аппроксимации основных уравнений. Распределение давления находится численным решением уравнения Пуассона. Выполнены параметрические исследования кинематических и динамических характеристик течения в зависимости от основных параметров задачи для неньютоновских сред. Продемонстрировано влияние числа Рейнольдса, показателя нелинейности реологической модели и геометрических характеристик препятствия на значение коэффициента местного гидравлического сопротивления.
Традиционно трубы с препятствиями в виде локальных участков сужения / расширения применяются в качестве конструктивных элементов современного теплообменного оборудования. Кроме того, исследование течений жидкости в таких каналах проводится в связи с изучением сосудистых патологий в кровеносных системах.
Для оценки влияния препятствия на распределения кинематических и динамических характеристик течения проведены многочисленные экспериментальные и численные исследования потоков в каналах и трубах с различной формой препятствия. Геометрия препятствия, рассматриваемая в данной работе, часто применяется при моделировании течений крови в стенозированных сосудах. Исследованию течений крови в сосудах с использованием ньютоновской модели реологического поведения посвящено большое количество работ [1–15]. Авторы экспериментальной работы [1] изучили течение жидкости в трубе для широкого диапазона чисел Рейнольдса и геометрических параметров препятствия. По результатам этих экспериментов удалось установить критические значения числа Рейнольдса для заданных геометрий, при которых происходит переход режима течения от ламинарного к турбулентному, описать структуру потока и определить потери давления, связанные с наличием препятствия в трубе. Известны исследования, например [4–7], где численно моделируется стационарный поток крови через стенозированную артерию. В [4] форма стеноза задается косинусоидальной функцией, а в [5, 6] – гауссовой функцией. Работа [8] посвящена изучению влияния формы стеноза на структуру течения. Авторы рассматривают четыре формы стеноза: прямоугольную, трапецеидальную и стеноз, задаваемый функцией косинуса и Гаусса. В [9] исследуются течения в каналах с несимметричными конфигурациями препятствия. В работах [10, 11] моделируются течения жидкости через два последовательно расположенных сужения, изучено влияние наличия второго сужения на распределения характеристик потока. В [12–14] рассматривается нестационарное течение ньютоновской жидкости в канале с препятствием заданной формы.
Экспериментальные исследования показали, что кровь при низких скоростях сдвига проявляет неньютоновские свойства. Имеется значительно меньше работ по моделированию течений в каналах с препятствием, где в качестве исследуемой жидкости рассматривается реологически сложная среда [16–21]. В [16] численно исследуется течение в крупной артерии, а в качестве рассматриваемой среды используется вязкоупругая жидкость. В [17] для описания реологических свойств течения жидкости в плоском канале авторы используют степенную модель. В [18–20] исследуются пульсирующие течения неньютоновской среды, для описания которой применяются различные реологические модели, в частности, в [18] используется модель Балкли–Гершеля. В отмеченных работах по численному моделированию течений используется двумерная постановка задачи. Развитие средств численного моделирования позволило авторам [22] рассмотреть трехмерное течение в круглой трубе с препятствием, используя пакет прикладных программ ANSYS.
Цель настоящей работы заключается в изучении характеристик потока степенной жидкости в трубе с препятствием заданной формы и в расчете коэффициента местного гидравлического сопротивления для умеренных чисел Рейнольдса при различных значениях показателя нелинейности и геометрических характеристик препятствия.
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} }}$Реологические свойства среды описываются уравнением Освальда-де Ваале [24], в котором безразмерная эффективная вязкость вычисляется по формуле
Здесь 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 реализуются мягкие граничные условия. Таким образом, условия на границах в переменных функция тока – вихрь запишутся в виде
Граница Г2 описывается функцией
Значения 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 }$Область решения покрывается квадратной сеткой с шагом h. Для разностного представления уравнений в преобразованной системе координат (2.1) и (2.2) применяется явная разностная схема, конвективные слагаемые в основных уравнениях аппроксимируются разностями против потока. Условия сходимости итерационного процесса записываются в виде
Для подтверждения аппроксимационной сходимости выполнены расчеты на последовательности вложенных сеток. Рисунок 2 демонстрирует профили аксиальной скорости на последовательности сеток в сечении 2 (рис. 1). В табл. 1 приведены значения аксиальной скорости на оси симметрии в сечении Г3 для разных значений h. Результаты расчетов, представленные на рис. 2 и в табл. 1, демонстрируют аппроксимационную сходимость предложенного алгоритма. Все дальнейшие расчеты проводились на сетке с шагом h = 0.025.
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} $Для записи граничного условия на границе Г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}$.
Уравнение Пуассона для давления решается методом установления. Все пространственные производные аппроксимируются со вторым порядком точности. Для определения сходимости процесса используется неравенство
где ε принимается равным 10–7.4. РАСЧЕТ МЕСТНОГО ГИДРАВЛИЧЕСКОГО СОПРОТИВЛЕНИЯ
Перепад давления между входным и выходным сечениями трубы с препятствием можно представить в виде суммы [26]
Таким образом, коэффициент местного гидравлического сопротивления рассчитывается по формуле
Значение ∆pм вычисляется экстраполяцией давления с участков одномерного течения перед и за препятствием в сечение 1 (рис. 1). Давление на этих участках изменяется по линейному закону.
5. РЕЗУЛЬТАТЫ РАСЧЕТОВ
Согласно результатам [21], картины стационарного течения степенной жидкости в трубе с участком сужения/расширения характеризуются зонами одномерного течения перед и за препятствием и зоной двумерного течения в окрестности препятствия с образованием циркуляционной зоны. На рис. 3 показаны типичные распределения линий тока для псевдопластичной жидкости (m < 1), демонстрирующие появление циркуляционной зоны за препятствием с увеличением числа Рейнольдса.
На рис. 4а представлены зависимости перепада давления от числа Рейнольдса для ньютоновской жидкости. Такая область решения использовалась в работах [10, 27], где представлены результаты численного моделирования. Экспериментальное исследование течения в трубе данной геометрии описано в [1], где установлено, что значение числа Рейнольдса выше 250 при заданных геометрических параметрах приводит к потере устойчивости ламинарного режима. Наблюдается хорошее согласование результатов расчетов с данными, представленными в [1, 10, 27].
На рис. 4б показаны изменения давления вдоль оси симметрии трубы для трех значений показателя нелинейности. Полученное распределение давления при m = 1 практически совпадает с данными из [11]. Максимальное значение перепада давления реализуется в случае дилатантной жидкости (m = 1.2). Анализ графиков, представленных на рис. 4б, позволяет сделать вывод, что для псевдопластичной жидкости (m = 0.8) перепад давления между входным и выходным сечениями трубы уменьшается по сравнению с течением ньютоновской жидкости. Такое поведение распределения давления согласуется с изменениями распределений эффективной вязкости в зависимости от параметра нелинейности, представленными в [21].
На рис. 5 показаны распределения давления для псевдопластичной жидкости при трех значениях параметра α. Перепад давления в областях одномерного течения постоянный и соответствует течению Пуазейля в круглой трубе. Увеличение глубины перекрытия способствует росту зоны двумерного течения за препятствием. Формирование циркуляционной зоны за препятствием [21] приводит к возникновению области с минимальным давлением, локализованной в окрестности максимального сужения трубы.
Влияние параметра нелинейности для двух значений Re на коэффициент местного гидравлического сопротивления представлено на рис. 6. С уменьшением m коэффициент местного гидравлического сопротивления падает, что обусловлено ростом среднего уровня вязкости в области двумерного течения, влияние показателя нелинейности усиливается с уменьшением числа Рейнольдса.
Зависимости См от числа Рейнольдса для разных значений m и геометрических параметров α, Z0 представлены на рис. 7. Графические зависимости показывают, что с ростом числа Рейнольдса коэффициент местного гидравлического сопротивления монотонно уменьшается. Уменьшение глубины перекрытия канала α также приводит к уменьшению значений коэффициента местного гидравлического сопротивления. Влияние параметра Z0 становится существенным при малых значениях α.
ЗАКЛЮЧЕНИЕ
Выполнено математическое моделирование течения степенной жидкости в трубе с препятствием. Разработан и реализован численный алгоритм решения поставленной задачи. Проведены параметрические исследования динамических характеристик потока для показателя нелинейности $0.6 \leqslant m \leqslant 1.4$ и числа Рейнольдса в диапазоне от 1 до 100. Вычислены значения коэффициента местного гидравлического сопротивления и выявлены тенденции его изменения при исследовании основных параметров задачи. Установлено, что с увеличением числа Рейнольдса коэффициент сопротивления уменьшается. Рост показателя нелинейности способствует увеличению коэффициента cопротивления. Уменьшение глубины перекрытия канала приводит к уменьшению значения коэффициента.
Исследование выполнено за счет гранта РНФ (проект № 18-19-00021).
Список литературы
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.
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.
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&eid=2-s2.0-0027093633
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.
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.
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.
Pontrelli G. Blood flow through an axisymmetric stenosis // P. I. Mech. Eng. H. 2001. V. 215. № 1. P. 1–10.
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.
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.
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.
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.
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.
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.
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.
Егоров В.А., Регирер С.А., Шадрина Н.Х. Особенности пульсирующего течения крови через резистивные кровеносные сосуды // Изв. РАН. МЖГ. 1994. № 2. С. 83.
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.
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.
Tu C., Deville M. Pulsatile flow of non-Newtonian fluids through arterial stenosis // J. Biomech. 1996. V. 29. № 7. P. 899–908.
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.
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.
Рыльцев И.А., Рыльцева К.Е., Шрагер Г.Р. Кинематика течения степенной жидкости в трубе переменного сечения // Вестн. Том. гос. ун-та. Математика и механика. 2020. № 63. С. 125–138.
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.
Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 616 с.
Ostwald W. Ueber die rechnerische Darstellung des Strukturgebietes der Viskosität // Kolloid Zeitschrift. 1929. V. 47. № 2. P. 176–187.
Годунов С.К., Рябенький В.С. Введение в теорию разностных схем. М.: Физматгиз, 1962. 340 с.
Идельчик И.Е. Справочник по гидравлическим сопротивлениям / под ред. М.О. Штейнберга. 3-е изд. М.: Машиностроение, 1992. 672 с.
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.
Дополнительные материалы отсутствуют.
Инструменты
Известия РАН. Механика жидкости и газа