Numerical simulation of a power-law fluid flow in a channel with double constriction.pdf Поле течения неньютоновской жидкости в цилиндрических каналах в окрестности препятствий представляет интерес для исследователей, изучающих механику жидкостей. В инженерных приложениях подобные каналы используются в качестве комплектующих элементов различного рода теплообменных установок и гидравлических систем. Кроме того, в биомеханике течение вязкой жидкости в трубке с участком сужения / расширения, описываемым, например, функцией косинуса или Гаусса, применяются для моделирования течения крови в стенозиро-ванном сосуде. Эти обстоятельства обуславливают актуальность экспериментального и численного исследований течений в каналах с препятствиями заданной конфигурации. В настоящее время доступно большое количество численных и экспериментальных работ, в которых исследуются течения ньютоновской жидкости в канале с одним препятствием [1-6]. Результаты этих работ содержат сведения о структуре потока и об основных параметрах задачи, влияющих на характеристики течения. Например, в работе [2] авторы одними из первых численно исследовали влияние числа Рейнольдса и геометрических параметров на распределения основных характеристик задачи. Развитие вычислительных технологий позволило изучить течения неньютоновских жидкостей [7-12]. Авторы работ [7, 8] используют реологические законы Балкли - Гершеля, Каро, Оствальда - де Ваале, Кассона для оценки влияния неньютоновских свойств при моделировании течений в артериях. В обзорной части [13] приведены работы авторов, исследующих ламинарные стационарные и нестационарные течения, как для ньютоновской, так и неньютоновской жидкости в случае, когда в канал включено одно препятствие. 1 Исследование выполнено за счет гранта Российского научного фонда (проект № 18-19-00021). Численное моделирование течения степенной жидкости в канале с двойным сужением 77 Течению в трубе с несколькими последовательно размещенными участками сужения / расширения уделено меньше внимания, при этом во всех исследованиях в качестве рассматриваемой среды используется ньютоновская жидкость. Подобные трубки встречаются в системах охлаждения летательных аппаратов и ускорителей ракет, которые функционируют при вхождении в плотные слои атмосферы, где рабочим телом служит газ или жидкость. В [14, 15] изучено влияние числа Рейнольдса в диапазоне от 5 до 200 на распределения напряжения и давления, сделаны выводы о локализации максимальных значений вихря на стенке. Представлены картины течения, демонстрирующие поля вихря и функции тока в зависимости от основных параметров задачи, дана оценка влиянию второго препятствия на структуру потока и на параметры течения. В статье тех же авторов [16] расширен диапазон исследуемых чисел Рейнольдса до 400. В [17] описывается трехмерное моделирование стационарного течения жидкости через несколько последовательно расположенных препятствий. Решение задачи осуществляется с помощью модифицированного метода LBGK, в основе которого лежит метод решеточных уравнений Больцмана (LBM), где столкновение частиц учитывается моделью Батнагара - Гросса - Крука (BGK). Авторы ограничились исследованиями в диапазоне чисел Рейнольдса от 10 до 150. К основным выводам данной работы можно отнести следующее: метод LBGK является полезным инструментом для моделирования установившихся течений; результаты, полученные с помощью трехмерной модели, хорошо согласуются с данными из работы [14], где применяется двумерная постановка задачи. Целью данной работы является исследование структуры потока и характеристик течения неньютоновской жидкости в осесимметричном канале с двумя препятствиями в зависимости от параметров задачи и оценка влияния геометрических параметров одного из препятствий на распределения характеристик потока около другого. Постановка задачи Рассматривается стационарное течение степенной несжимаемой жидкости в круглой трубе с двумя последовательно размещенными в ней участками сужения / расширения, форма которых описывается функцией косинуса. Область течения представлена на рис. 1. Рис. 1. Область течения Fig. 1. Flow area Для математического описания процесса используется постановка задачи в переменных вихрь - функция тока [18]. Реологическое поведение жидкости описы- И.А. Рыльцев, О.Ю. Фролов, Г.Р. Шрагер 78 вается законом Освальда - де Ваале [19]. Система уравнений имеет вид Л dz2 д(мю) д(ѵю) _ д дг dr Re д2ю д2ю 1 дю ю df2 r dr Re д 2д d дѵ 5м д + 2 дд дю + 2 дд дю + dzdr I дг дz ) дz дz дг дг ^д2д д2д Ѵ дм дѵд ю дд дz2 дг2 А дг дг J r дг (1) () д2 ш д2ш 1 дш дr 2 r дr -г- + -21---- _ -гю ; дг 2 д _ A m-1 (3) Безразмерные компоненты скорости и вихрь в уравнении (1) определяются формулами 1 дш 1 дш дѵ дм r дr r дг дг дr Здесь м, ѵ - аксиальная и радиальная компоненты скорости соответственно, Re _ к - число Рейнольдса, р - плотность, U - среднерасходная скорость, A _ (2еі]-е]і )2, ву - компоненты тензора скоростей деформаций, к - консистенция жидкой среды, m - степень нелинейности жидкости. В качестве масштабов обез-размеривания приняты следующие величины: длины - радиус трубы r0; скорости - среднерасходная скорость U. Граничные условия в переменных вихрь - функция тока имеют вид [18] Г, : ш _ f urdr, ю _ --, при г _ 0; о д „ 1 д2ш .... ч Г2 : ш _ const, ю _---, при r _ f (г); r дп дш „ дю Г3 :-L _ 0, - _ 0, при г _ L; дг дг Г4 : ш _ 0, ю _ 0, при r _ 0, дм (4) где п - нормаль к границе Г2, м _ 3т +1 т +1 Граница Г2 задается функцией f (г) _ 1 - rm 1+1 д 1 - SilL I 1 - cos П( г Ll) а1 2 а2 2 г e[Ll; L2 ]; 1 “'2 I і П(г L3) 1 -21 1 - cos-- г е [3 ; L4 1, г е [0;L1)u (L2;L3)и (L4;L]. Численное моделирование течения степенной жидкости в канале с двойным сужением 79 Метод решения Для получения численного решения сформулированной задачи применяется метод установления, который позволяет получить стационарное решение. В этом случае организуется нестационарный процесс, решение которого с течением времени оказывается независимым от него и устанавливается к решению исходной стационарной задачи [20]. Для реализации нестационарного процесса в уравнения для вихря (1) и функции тока (2) добавляются производные по времени функций ю, у соответственно. Итерационный процесс продолжается до установления в пределах заданной точности. Для аппроксимации дифференциальных уравнений область течения с криволинейной границей f(z) трансформируется в прямоугольную с помощью введения новой системы координат £ = z, п = r / f (z). В координатах £, п уравнение переноса вихря (1) и уравнение Пуассона для функции тока (2) с введенными производными по времени принимают вид дю д (ию) д(ию) 1 д(ѵю) dt д£ ^ дп f дп где Re ^ д2ю тт дю , -- + H--2 gд£2 дп д£дп д2 ю + G д 2ю 1 дю ю дп2 f2П дП f2П2 Re +£.*+|н -_! дt д£2 V f п) дп ду т д2у + „ д2у - 2 g-- + G- = - /пю, д£дп дп2 (5) (6) ,г_ - г_ ' f2 f ; g = п f_ f' S = 2 S = S1 + S2 + S3 + S4, ( Л д2.. „ з.. „ V 1 дѵ ди 1 д2ц g дц g д2ц f д£дп /п дп f дп2 ди f дп д£ g дп s2 = 2 - g ^ - g ^ 1 + 2± ^ ^, ч д£ дп J V д£ дп ) f2 дп дп S3 = Ѵц.+н ф - 2 g д£2 дп д£дп дѵ S4 = f2 )дп2 ю дц 1 ди дѵ 7"дп + "7 g дп f2 п дп Численное решение дифференциальных уравнений осуществляется конечноразностным методом. Преобразованная прямоугольная область решения покрывается равномерной в каждом направлении разностной сеткой пн ={£г = 'К Ч} = jh, i = 0,...,^ j = 0,...,N2}, где h - шаг сетки; N - количество узлов в направлении £; N2 - количество узлов в направлении п. Разностное представление уравнений (5), (6) выполняется с использованием явной разностной схемы. Конвективные слагаемые в уравнении (5) аппроксимируются схемой против потока. Условие сходимости вычислительного процесса к стационарному решению задачи имеет вид И.А. Рыльцев, О.Ю. Фролов, Г.Р. Шрагер 80 max ?. І < E. max ?. i 1 - < < E. где e - параметр сходимости; q - номер шага по фиктивному времени. Значение параметра e определяется в ходе численного эксперимента (e = 10-5). Выбор явной схемы обусловлен простотой реализации расчётного алгоритма. Схема против потока для аппроксимации конвективных слагаемых обеспечивает устойчивость разностной схемы. При расчете эффективной вязкости в случае. когда m < 1 значения ц стремятся к «бесконечным». Для устранения этой особенности используется модифицированная запись реологического уравнения (3). которая имеет вид ц = (A + X )m-1. (7) Здесь X - параметр регуляризации. значение которого выбиралось экспериментально и принималось равным 0.001. Для тестирования вычислительного алгоритма выполнены расчеты на вложенных сетках. В таблице приведены значения аксиальной скорости на оси симметрии в трех контрольных сечениях в зависимости от шага сетки. демонстрирующие аппроксимационную сходимость алгоритма. Все дальнейшие расчеты проводились с использованием шага 0.025. Аппроксимационная сходимость при Re = 10, m = 0.9, a1 = a2 = 0.5, Z1 = Z2 = 1, L1 = 3, L3 = 7 h umax. Z = L]+Z] umax. z=L3+Z2 umax. z = L 0.1 5.3844 5.4860 1.9408 0.05 5.4081 5.49675 1.9469 0.025 5.415 5.4956 1.9478 0.0125 5.4159 5.4959 1.9476 На рис. 2 показано сравнение распределений функции тока. полученных в настоящей работе. с данными. взятыми из работы [14]. при этом функция f(z) задаётся в соответствии с [14]. Структура потока характеризуется зонами одномерного течения в окрестности входной и выходной границ. В области препятствий течение имеет двумерный характер. за препятствиями образуются циркуляционные зоны. Наблюдается качественное согласование результатов. а т-1-1-1-1-1- 2 4 6 8 10 12 Рис. 2. Сравнение распределений линий тока из статьи [14] (а) и настоящей работы (b) (Re = 25. m = 1. a! = a2 = 0.5. Z\\ = Z2 = 1) Fig. 2. Comparison of streamline distributions (a) from the paper [14] and (b) in the present work (Re = 25. m = 1. a! = a2 = 0.5. Z\\ = Z2 = 1) Численное моделирование течения степенной жидкости в канале с двойным сужением 81 Восстановление давления Постановка задачи в переменных вихрь - функция тока не содержит давления. Для расчета динамических характеристик потока воспользуемся уравнением Пуассона [18], которое позволяет по известным полям скорости, вихря и функции тока рассчитать поле давления и в физических переменных имеет вид Re 92 p 92 p 1 dp -U +-Г- +-- = 2 9z2 dr2 r dr v 9ц 2 9ц 2 92ц - ■V2u + -■ V v + ■ ^ du dv dz dr du dv dr dz r 2 2 du dv J d 2ц du d 2ц dv dz dr dzdr V dr dz J dz2 dz dr 2 dr -a \\ (8) Для обезразмеривания величиныp используется масштаб pU. Уравнение (8) в преобразованной системе координат (t, ц) принимает вид Г Л \\ Г / \\2 Л Гд2 p d2p 2 g-Р- + dt2 д^дц H + - 1 dp+G ddu дц дц2 = 2 S5 ■ S6 -S7 ■ S8-f Re S9 V S10 + S11 + fj + S12 V S13 + S14 + уц^) + S15 + S7 ) + S16 S5 + S17 S6 ) ; (9) S5 = 9u 9u 0 1 9v 0 9v дv -g-, S6 =--, S7 =--g- 9t дц f дц 9t дц 92u „ д2 u S10 =W - 2g „ 1 9u дц дц S8 = , S9 = --g -, f дц dt дц +Hdu + 2 ddu s = -d-Us = -дц д^дц дц g дц2, 11 f2 дц2 , 12 f дц, д 2 v S'3=dF- 2 g d2v . dv д2 v 1 d2v 1 д2ц g дц g д2ц + H--+ g -~, S14 = - ---, S15 = д^дц дц дц2 f2 дц2 f д^дц f\\ дц f 9ц2 S16 =92т - 2 g Ц + H £ + g2- 9±. -2 д^9ц 9ц ^2 dV f2 дц2. S17 = 9t2 д^дц дц дц2 Для постановки граничного условия на границе Г2 используются уравнения движения, которые с учетом условий прилипания на твердой стенке при ц = 1 в системе координат (t, ц) запишутся в виде 9p дц f 9p 1 Г ц Гдю ю | ю дц |; дц Re Vf 1дц+ц ) + f дц ); (10) Г Гдю дю Л Г дц дц Л Л 5Іц Vsf - g Зц)-ю Ы- g st )J. (11) 9p df Используя (10) и (11), получим уравнения для расчета давления на границах Г1 и Г2. В выходном сечении Г3 безразмерное давление задается равным нулю, на оси симметрии Г 4 реализуется условие симметрии. Таким образом, граничные условия для расчета давления принимают вид Г • 9? = 1 Г ц Г дю + ю| + ю дц 1: dt=-Re V7 V^+f)+f сц Г 9p 1 Г J Гдю дю dt=Re Igf Hdf- g de 2 - ю при t = 0; дц дц g ц Г дю + юЛ + ю дц 9t дц)) V f Ѵдц f) f дц приц = 1; И.А. Рыльцев, О.Ю. Фролов, Г.Р. Шрагер 82 Г3:p = 0 при^ = L; dp Г4:- = 0 приц = 0. дц Решение дифференциального уравнения для давления, как и в случае нахождения функций ю, у, осуществляется конечно-разностным методом в явном виде с использованием метода установления. Условие для определения сходимости по давлению имеет вид max и і 1 - РІ, і < е, где е = 10 7. Результаты Рисунок 3 демонстрирует типичные картины течения псевдопластичной жидкости (т < 1) в виде распределений линий тока в зависимости от Re. При числе Рейнольдса равном 5 за препятствиями формируются циркуляционные зоны одинакового размера. При дальнейшем увеличении Re до 15 циркуляционная зона за первым участком сужения / расширения увеличивается и занимает всю область между препятствиями выше сечения r = 0.5, зона за вторым препятствием также увеличивается. На рис. 4 второе препятствие смещено относительного первого в сторону выходной границы, в результате тенденция распределения линий тока для Re, равного 1 и 5, сохраняется, а для Re = 15 зона за первым перекрытием стала равной зоне за вторым. Рис. 3. Линии тока (т = 0.8, а! = а2 = 0.5, Z1 = Z2 = 1) Fig. 3. Streamline distributions (m = 0.8, а! = а2 = 0.5, Z\\ = Z2 = 1) Рис. 4. Линии тока (m = 0.8, а! = а2 = 0.5, Z\\ = Z2 = 1) Fig. 4. Streamline distributions (m = 0.8, а! = а2 = 0.5, Zi = Z2 = 1) Численное моделирование течения степенной жидкости в канале с двойным сужением 83 Влияние параметра нелинейности на структуру потока показано на рис. 5. С уменьшением m размер циркуляционных зон за препятствиями увеличивается, начиная с некоторого значения параметра нелинейности, первая зона занимает все пространство между препятствиями, а размер второй продолжает расти в направлении выходной границы. Рис. 5. Линии тока (Re = 25, а! = а2 = 0.5, Z\\ = Z2 = 1) Fig. 5. Streamline distributions (Re = 25, а1 = а2 = 0.5, Z1 = Z2 = 1) На рис. 6 показано влияние глубины второго перекрытия (а2) на распределение линий тока в случае течения псевдопластичной жидкости. Уменьшение а2 приводит к распространению первой циркуляционной зоны в область второго препятствия. Рис. 6. Линии тока (Re = 25, m = 0.8, а1 = 0.5, Z1 = Z2 = 1) Fig. 6. Streamline distributions (Re = 25, m = 0.8, а1 = 0.5, Z1 = Z2 = 1) Графики на рис. 7 демонстрируют сравнение распределений давления на твердой стенке с данными [16]. Для сравнения динамических характеристик потока функция f(z) задавалась в соответствии с [16], кроме того, безразмерное давление на входе в канал принималось равным нулю. В зоне одномерного течения давление падает линейным образом, а на участках с препятствиями распределение давления носит сложный характер. Кривые на рис. 7 демонстрируют согласование полученных результатов с данными [16]. На рис. 8 показаны распределения давления вдоль оси симметрии для трех значений параметра m в канале, изображенном на рис. 5. Минимальному перепаду давления между входом и выходом соответствует течение вязкопластичной жидкости, что согласуется с изменением эффективной вязкости в зависимости от m. И.А. Рыльцев, О.Ю. Фролов, Г.Р. Шрагер 84 Рис. 7. Распределения давления вдоль стенки (Re = 2.5, m = 1, aj = a2 = 0.333, Zj = Z2 = 1). Обозначения:--настоящая работа; •-•-• - работа [16] Fig. 7. Pressure distributions along the wall (Re = 2.5, m = 1, a1 = a2 = 0.333, Z1 = Z2 = 1). Notations:--present work; •-•-• - paper [16] Рис. 8. Распределение давления вдоль оси симметрии (Re = 25, a! = a2 = 0.5, Z! = Z2 = 1). Обозначения: - m = 0.8,--m = 1,----m = 1.2 Fig. 8. Pressure distributions along the symmetry axis (Re = 25, a! = a2 = 0.5, Z! = Z2 = 1). Notations:.....- m = 0.8,--m = 1,----m = 1.2 На рис. 9 и 10 приведены распределения аксиальной скорости на оси симметрии канала для псевдопластичной, ньютоновской и дилатантной жидкостей. Графики на рис. 9 и 10 показывают, что рост числа Рейнольдса приводит к уменьшению максимальных значений скорости в зоне двумерного течения, что согласуется с кинематикой течения, показанной на рис. 3 и 5. Рисунок 11 иллюстрирует распределения вихря на твердой стенке в зависимости от местоположения второго препятствия относительного первого для трех значений параметра Re. Расчеты показали, что увеличение числа Рейнольдса приводит к росту максимальных значений вихря на стенке. С изменением геометрии происходит соответствующее перераспределение максимумов в функции ra(z). 0 4 8 12 z 0 4 8 12 z 0 4 8 12 16 20 z Рис. 9. Распределения аксиальной скорости на оси симметрии (а1 = а2 = 0.5, Z1 = Z2 = 1, L1 = 3, L3 = 5). Обозначения: - m = 0.8,--m = 1,----m = 1.2 Fig. 9. Axial velocity distributions along the symmetry axis (a1 = a2 = 0.5, Z1 = Z2 = 1, L1 = 3, L3 = 5). Notations:.....- m = 0.8,--m = 1,----m = 1.2 u Рис. 10. Распределения аксиальной скорости на оси симметрии (а1 = а2 = 0.5, Z1 = Z2 = 1, L1 = 3, L3 = 7). Обозначения:.....- m = 0.8,--m = 1,----m = 1.2 Fig. 10. Axial velocity distributions along the symmetry axis (a1 = a2 = 0.5, Z1 = Z2 = 1, L1 = 3, L3 = 7). Notations:.....- m = 0.8,--m = 1,----m = 1.2 Рис. 11. Распределения вихря на стенке (m = 1, Re = 5, а1 = а2 = 0.5, Z1 = Z2 = 1). Обозначения: - L1 = 3, L3 = 5;----L1 = 3, L3 = 7;--L1 = 3, L3 = 11 Fig. 11. Shear stress distributions on the wall (m = 1, Re = 5, a1 = a2 = 0.5, Z1 = Z2 = 1). Notations: - L1 = 3, L3 = 5;----L1 = 3, L3 = 7;--L1 = 3, L3 = 11 И.А. Рыльцев, О.Ю. Фролов, Г.Р. Шрагер 86 Заключение Выполнено численное моделирование течения неньютоновской жидкости в осесимметричном канале с двумя препятствиями. Создана и протестирована программа расчета течений для ЭВМ. Показаны и описаны картины течения степенной жидкости в зависимости от параметров процесса. Проведены параметрические расчеты динамических и кинематических характеристик течения.

Скачать электронную версию публикации