CC BY
19
УДК 532.546+532.55
DOI:10.30724/1998-9903-2019-21-5-97-109
КОНЦЕПЦИЯ ТУРБУЛЕНТНОЙ «ВИХРЕВОИ ЗАСЫПКИ» - МОДЕЛИ И
МЕТОДЫ
Л.Э. Меламед1, Г.А. Филиппов'
2
Акционерное общество «Интеллект», г. Москва 2Отделение энергетики, машиностроения, механики и процессов управления Российской академии наук, г. Москва
melamed@yandex. гы
Резюме: Представлены модели и методы изучения турбулентности, основанные на концепции турбулентной «вихревой засыпки». Суть этой концепции состоит в том, что турбулентное течение рассматривается как ламинарное, текущее через «вихревую засыпку», создающую внутреннее сопротивление. Это сопротивление можно рассматривать либо как распределенное, либо как локально-сосредоточенное. На основе первого представления получено модифицированное уравнение Навье-Стокса, его приближенное аналитическое и численные решения. На основе второго представления и разработанного для этих целей метода локальных флуктуаций получена компьютерная модель турбулентного потока в трубах. С помощью моделирования показано, что при задании определенной системы локальных флуктуаций вязкости расчетный профиль течения соответствует профилю скорости турбулентного потока. Величина и профиль турбулентной вязкости потока полностью определяются структурой и свойствами «вихревой засыпки». Результаты работы подтверждают возможность и эффективность рассмотрения турбулентности на основе данной концепции.
Ключевые слова: турбулентность, вихревая засыпка, моделирование, локальные флуктуации.
Благодарности: Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (РФФИ, проект № 18-08-00051а).
Для цитирования: Меламед Л.Э., Филиппов Г.А. Концепция турбулентной «вихревой засыпки» - модели и методы // Известия высших учебных заведений. ПРОБЛЕМЫ
ЭНЕРГЕТИКИ. 2019. Т. 21. №5. С. 97-109. doi:10.30724/1998-9903-2019-21-5-97-109.
THE CONCEPT OF TURBULENT «VORTEX BACKFILL» - MODELS AND
METHODS
1 Joint-stock company "Intelligence", Moscow, Russian Federation 2 Department of power, mechanical engineering, mechanics and control processes of Russian Academy of Sciences, Moscow, Russian Federation
lev. melamed@yandex. ru
Abstract: Models and methods for studying turbulence based on the concept of turbulent "vortex backfill" are presented. The essence of this concept is that the turbulent flow is considered as laminar, flowing through a "vortex backfill ", which creates internal resistance. This resistance can be considered either as distributed, or as locally concentrated. Based on the first representation, a modified Navier-Stokes equation, its approximate analytical and numerical solutions are obtained. Based on the second concept and the local fluctuation method developed for these purposes, a computer model of the turbulent flow in the pipes is obtained. Using simulation, it is shown that, when a certain system of local viscosity fluctuations is specified, the calculated flow profile corresponds to the profile of the turbulent flow velocity. The magnitude and profile of the turbulent viscosity of the flow are completely determined by the structure and properties of the "vortex backfill ". The results of the work confirm the possibility and efficiency of considering turbulence based on this concept.
LE Melamed1, GA Filippov2
2
Keywords: turbulence, «vortex backfill», simulation, local fluctuations.
Acknowledgments: This work is fulfilled with financial support of the Russian fund of basic researches (RFFI, the project № 18-08-00051a).
For citation: Melamed LE, Filippov GA. The concept of turbulent «vortex backfill» - models and methods. Power engineering: research, equipment, technology. 2019;21(5):97-109. (In Russ). doi:10.30724/1998-9903-2019-21-5-97-109.
Введение
Турбулентные течения жидкостей и газов присутствуют в огромном числе машин и аппаратов не только в энергетике (в частности, в атомной энергетике), но ив других отраслях техники и явлениях природы. Изучение турбулентности наталкивается на серьезные экспериментальные и математические трудности, но «... это не должно затемнять того факта, что сердцевиной процесса является физика» [1]. В настоящее время экспериментальное изучение турбулентного течения в трубах ведется с огромным размахом. Кроме специально созданных экспериментальных установок в Делфте, Голландия, Принстоне, США (установка Superpipe), были созданы огромные экспериментальные установки CoLaPipe в Коттбусе, ФРГ и CICLoPE в Форли, Италия.
Этому изучению способствуют и модели, физические и математические. Модели [2-5] заостряют внимание на основных аспектах физики процесса, оставляя в стороне второстепенные. В физике турбулентности основную роль играют вихри - вращающиеся объёмы вещества, образующиеся в потоке после достижения им некоторой критической скорости. Рассмотрению и учету вихревых структур посвящено огромное количество работ. В качестве одного из последних примеров можно привести работу [2], в которой вихри (или диполи) считаются квазитвердыми и обтекаются абсолютно невязкой жидкостью. В данной работе рассматривается другая модель, а именно, предложенная в [3] модель турбулентной «вихревой засыпки», в которой вязкость несущего потока учитывается в полном объёме. Концепция этой модели состоит в том, что турбулентный поток можно рассматривать как ламинарный, но текущий через «вихревую засыпку», которую нужно учитывать как внутреннее сопротивление, распределенное или локально-сосредоточенное. При этом используется глубокая внутренняя связь между турбулентными течениями и течениями в зернистом слое. На основе этой концепции предложены методы ее реализации и получены новые результаты.
Модель турбулентной «вихревой засыпки» построена как двухфазная. Течение в ней можно представить себе как ламинарный поток между жидкими вихрями, движущимися и вращающимися, меняющими форму и размеры. Поступательно они двигаются с чуть меньшей скоростью, чем свободный поток. Вращение вихрей делает их мало «прозрачными» для основного, несущего потока, что увеличивает их эффективную вязкость. Эти обстоятельства создают внутреннее сопротивление для потока, которое отсутствует в ламинарном течении.
Термин "засыпка" применяется к зернистым слоям из твердых частиц, применяющимся в атомной энергетике, химической и других отраслях промышленности [6]. В аппаратах с такими слоями потоки рассматриваются как преодолевающие внутренние объемные сопротивления. Величины таких сопротивлений находятся на основе экспериментальных данных.
Модель с равномерно распределенным сопротивлением и ее математическая реализация
Расчетные модели турбулентности, применяющиеся в настоящее время, используют, в основном, понятие о турбулентной вязкости. Внутреннее сопротивление рассматривается редко. Примером является работа [7]. В ней сила сопротивления добавлена в правую часть уравнения Навье-Стокса в виде бесконечного степенного ряда по скоростям. Коэффициенты ряда определяются по экспериментальным данным. В нашей работе добавлено только одно слагаемое. Экспериментальные данные не требуются.
Уравнение движения в круглой цилиндрической трубе имеет вид (здесь F - распределенное сопротивление)
1 d C^W) = Р - F или уч (! _ F). (1)
r dr dr x pwa P
Здесь введены следующие обозначения: w - скорость, м/с, y = w/wa - относительная скорость, wa - средняя скорость потока, м/с, r - радиус, м, x = r/R - относительный радиус,
Я - радиус канала, м, ц - динамическая вязкость, Па с, Р = др / дz - градиент давления, Па/м, р - давление, Па, г - направление по оси потока, м, Е - градиент давления от распределенного сопротивления, Па/м.
Покажем, что при определенном выражении для функции Е/Р = у(х,у) это уравнение определяет профиль осредненной скорости турбулентного потока. Величины Р (полное сопротивление турбулентного потока) и Е (сопротивление «вихревой засыпки») в общем виде равны выражениям
р 2 с 2
Р = Ар ^2 и р = к р ^ , (2)
2Я 2 2
где р - плотность, кг/м3, ds - характерный диаметр зерна «засыпки», м. Нас интересует входящее в уравнение (1) отношение у = Е/Р. Для определения этого отношения была проведена серия вариантных расчетов, в которых менялись как число Re, так и гипотетически возможные параметры зерен «вихревой засыпки». Коэффициент сопротивления потока 4 в его зависимости от числа Яе = 2Я^ар/д известен по экспериментальным данным и их аппроксимациям. Коэффициент сопротивления засыпки рассматриваемый в технических приложениях [4], зависящий от пористости засыпки и диаметра ее «зерна», также известен. Результаты расчетов показали, что по порядку величин градиенты давления Р и Е близки. Оказалось, что величина у хорошо описывается формулой у = (у/у.)2 , где ус = у(0)(\-х2/2)°'5 - профиль, аппроксимирующий скорость турбулентного потока в его центральной части. Таким образом, величина у, которая представляет собой долю общих потерь, приходящуюся на засыпку, зависит только от радиуса и формы профиля скорости и не зависит от ее амплитуды.
Теперь модифицируем уравнение (1). Заменим величину у выражением (у/ус)2 и обозначим РЯ2/^м>а= О. Введем коэффициент 9 = 9(х) - характеристическую функцию, значение которой определяется зоной течения. Коэффициент 9 = 1 в турбулентном ядре и 9 = 0 в вязком подслое, вблизи стенки. В переходной зоне, в первом приближении, положим величину 9 = 0,5. Таким образом, радиусы точек перемены значений 9 зависят только от числа Re.
Уравнение движения примет вид
2
—' + — = в(1 - 0(х)^-У—«-) , 0 < X < X '=1,2,3. (3)
х у2(0)(1 - х2/2)
Уравнение неразрывности в нашем случае заменяется условием постоянства объемного расхода
Я 2 1
2п | ^(г)гйг = м>апЯ или 21 yxdx = 1.
0 0 Это условие определяет полный расход через все сечение, а при заданном расходе -среднюю скорость
Для отдельных зон течения условие для расхода приобретет вид
х1 х1
2 / —^¿х = 2 / фХйХ = Т ,=1 2 3 (4)
0 0
где Ji - относительный расход, а ф(х) - принятое за эталон распределение скоростей, а именно - распределение по формуле Рейхардта, приведенной к виду ф(х) = w/wa.
Итак, задача в окончательном виде содержит уравнение движения (3), уравнение неразрывности (4) и граничные условия: на левой границе (при х = 0) у,(0) = 0, на правой границе (при х = х,) у = у * .
Расчет поля скоростей может вестись как для всех зон, последовательно от стенки к оси, так и для каждой зоны отдельно. В этом случае граничные условия для них можно получить из эталонного решения. В пристенной зоне, при 9=0, решение соответствует известному параболическому закону.
Решение одномерного уравнения
Рассмотрим нахождение профиля скорости в ядре турбулентного потока. Интегральное условие (4) заменяется эквивалентным дифференциальным (с введением новой неизвестной и дополнительным краевым условием), и задача становится системой трех обыкновенных нелинейных дифференциальных уравнений первого порядка. Величина у(0) определяется итерационно. Ее начальная величина берется из известного диапазона значений 1,15 - 1,30. Распределение ф определяется формулой Рейхардта. Коэффициент 9=1 в ядре потока. Величина О, на основе первой из формул (2) и знака градиента, равна О = -4(Яе)Яе/8. Решение производится с помощью системы ЫайаЪ.
Полученные результаты (профили скоростей в турбулентном ядре потока) представлены на рис.1. Сравнение расчетных и эталонных скоростей показывает их близость. Скорости на осях как при Re = 4-103, так и при Re = 4-104, оказались практически одинаковыми, при Re = 4 105 их расхождение не превысило 1,4%.
Включение внутреннего сопротивления в двумерную задачу
Эта же задача, но уже для двумерной осесимметричной задачи и для всего поля течения (а не только для турбулентного ядра потока), была решена в [3] с помощью компьютерной системы Сот^'о1 МиШрку81с8. В этой системе есть возможность в уравнении течения задать дополнительное сопротивление ¥ аналитически, а именно, в виде
= 9^--(5)
дг у(0)2 (1 - х / 2) где у, м/с - скорость в направлении оси г трубы.
Уравнение неразрывности не изменяется. Таким образом, проверено утверждение о том, что эта методика дает возможность рассчитывать турбулентные течения по программам для ламинарных течений.
Приближенное аналитическое решение
Одновременно с экспериментальными исследованиями [8-10] ведется и теоретический анализ процесса, связанный с выбором «правильной» формы аппроксимации получаемых экспериментальных данных [11-13], удовлетворяющей при этом определенным теоретическим представлениям. Обсуждение этой темы велось ранее крупнейшими учеными и ведется сейчас (Г.И. Баренблатт [14,16], И.И. Вигдорович [15] и др.) и зачастую принимает форму публичной дискуссии. Сталкиваются два закона изменения средней скорости по радиусу канала (точнее - в ядре потока, вне промежуточного слоя и ламинарного подслоя) - логарифмический (универсальный, не зависящий от числа Re) и степенной (не универсальный). Авторы [14,16] пишут, что «универсальный логарифмический закон не соответствует экспериментальным данным». Они предлагают новый степенной закон, включающий число Re. Автор [15] возражает. Предлагаемое на основе данной концепции уравнение турбулентного движения, как будет показано ниже, имеет степенное решение.
Y
0.6
0.4
0.6
0.8
Х
Рис. 1. Расчетные (пунктирные) и эталонные (сплошные) скорости турбулентного ядра потока в трубе при различных числах Re (х - относительный радиус, у - относительная скорость). Профили соответствуют Re = 104, 105, 106 - сверху вниз (в левой части рисунка)
Найдем приближенное аналитическое решение уравнения (1) для турбулентного ядра. Оно, как известно, при больших Re занимает практически все поле течения [17]. Граничную точку ядра (хк, ук) определяем в соответствии с трехслойной схемой Кармана. Кроме выше названных, введем еще новое условие, а именно - условие в особой точке. Как отмечено в [8] по результатам экспериментов, так и по обобщающим более ранние эксперименты выражениям [17] имеет место следующее. На графиках кривых у = fxДe), определяющих скорость, существует точка (х*,у*), в которой все профили (для всех Re) пересекаются, и при этом у*=1. Эту особую точку, не зависящую от числа Re, мы впервые предлагаем использовать в качестве ещё одного условия задачи:
у(х*)=у*, (6)
Это условие позволит определить дополнительный параметр решения. Получим приближенное аналитическое решение, применив метод последовательных приближений. Как известно, успех этого метода зависит от подходящего выбора начального (нулевого) приближения. Это нулевое приближение примем в виде:
У = У(0).
1 п\ 1 - с - х I
(7)
Интегрирование уравнения (1) с граничными условиями на оси и при х = хк (при 9=1) приводит к степенному решению (первому приближению)
2
У = ук
■О
V х2 ) + (хк'хР)/Р
(8)
где т = с/4, р = п+2. Условия (4) и (6) позволяют найти неизвестные параметры т и р, т. е. полностью определить параметры скорости.
Приведем некоторые подробности вычислений. После подстановки решения (8) в (4) условие баланса расхода станет таким
р+2
Ук. 2
т
о-+-
2 Р( Р + 2)
<2ъ
- = 0
Условие в особой точке станет таким:
у* - Ук - О т (х2 - хк ) + (
Р Р - х ^
/ р
= о
Входящие в эти два уравнения параметры т и р определяем следующим образом. Коэффициент сопротивления примем равным £ = (1,821§Ке-1,64)"2 (известная формула Г.К. Филоненко). Особая точка по Рейхардту [17] имеет координаты х*~0,80, у*=1. Далее вспомним тождество О = -£, Яе/8, которое следует из определения коэффициента сопротивления £ с учетом знака градиента давления Р. Полученное решение достаточно близко совпадает с экспериментальными данными для ядра потока. Решение (8) удобно представить в виде
+ а (х1- х 2 ) + ь (хк- ^ )',
У = У к
(9)
где а = -От, Ь = -О/р . Вычисленные нами для диапазона 4 10 ^е<4 10 параметры а, Ь и р описываются следующими выражениями:
а = 0,0175^2Яе - 0,27ЩЯе + 1,285, Ь = 0,0841^е - 0,0865, р = 0,373ехр (0,8511^е).
Проанализируем качественные характеристики решения. Рассмотрим, из каких компонентов складывается модуль скорости и как они меняются в зависимости от числа Re. Этих компонентов три: ук - «несущая» скорость турбулентного ядра; У1=а(хк2-х2) -параболическая составляющая скорости (она создает профиль турбулентного ядра); У2=Ь(хкр-хр) - степенная непараболическая составляющая скорости, которая оказывает очень малое влияние на профиль потока всюду, за исключением пристеночной области. На рис.2 представлены эти компоненты при двух числах Яе - при Re=1•104 (рис.2а) и Re=4•107 (рис.26).
о.в 0.6 0.4 0.2
Ук
У1
у2
О.В 0.6 0.4 0.2
У2
Ук
Уг
0.5 х
а)
0.5 х
Ь)
Рис.2. Составляющие скоростей при Re=1■104 (а) и Re=4■107 (Ь); здесь ук - несущая скорость, у1 - параболическая компонента скорости, у2 - степенная (непараболическая) компонента скорости
Несмотря на изменение чисел Re более чем на три порядка, характер компонент сохраняется неизменным. Величина ук постоянна по радиусу; величина >>1 - парабола с переменными параметрами; величина >2 почти постоянна по радиусу и уменьшается только вблизи края зоны ядра или стенки; при этом все они зависят от Re. Коэффициент а в рассматриваемом диапазоне чисел Re уменьшается от 0,45 до 0,2, коэффициент Ъ растет от 0,25 до 0,55. Сильнее всего меняется показатель степени р - от 11 до 240. При таких больших степенях р значение хр близко к нулю почти всюду, что и приводит к почти полному постоянству >2 по радиусу.
Увеличение параметра Ъ и координаты хк с ростом Re приводит к увеличению компоненты >2 почти в шесть раз. Однако самое интересное происходит с компонентой у1. Она уменьшается мало, только в два раза, и ее вклад остается существенным при самых больших числах Рейнольдса. Именно ей обязан своей параболичностью профиль скорости в его большей части при любом значении числа Re. Итак, полученное решение и его анализ позволяют оценить форму и величины компонент, формирующих профиль турбулентной скорости.
Важно отметить, что выражение (9) с условиями (4) и (6) пригодно для описания точного решения для ламинарного течения. Тогда константы будут такими: хк=1, >к=0 (условия на стенке), Qк=1 (полный расход), >*=1, х*=1/^2 (особая точка в данном случае). Подставив эти величины в условия (4) и (6), получим систему, откуда найдем т=0, р=2. Для этих значений из (9) получим у = О(1-х2)/4. Так как для ламинарного потока 4=64/Яе, получим окончательно у= 2(1-х2), т.е. скорость ламинарного потока. Таким образом, формула (4) является обобщающей для турбулентного и ламинарного режимов течения.
В завершение этого раздела можно констатировать, что полученная нами степенная, зависящая от числа Яе форма профиля турбулентного течения в трубах представляется нам предпочтительной в дискуссии [14-16].
Модель с локально-сосредоточенными сопротивлениями и ее компьютерная реализация
Методы математического моделирования не дают возможности явным образом вводить в рассмотрение основные физические составляющие турбулентного течения, а именно, сами вихри или вихревые трубки. Их воздействие заменяется некоторыми распределенными внутренними параметрами (турбулентной вязкостью, кинетической энергией, скоростью диссипации кинетической энергии и т.п.). Явное же введение в компьютерные расчеты этих элементов крайне затруднительно (в связи с их многочисленностью, необходимостью задания граничных условий). Предложенный нами метод локальных флуктуаций [16] дает возможность обойти вышеназванные затруднения. Рассмотрим его.
Метод локальных флуктуаций
Сутью метода является замена вихрей эквивалентными им (по динамическому воздействию) флуктуациями вязкости. В рассмотрение вводится задаваемая априори система вихрей, которая рассчитывается впрямую. Работа является продолжением ряда работ авторов, в которых рассматривалось влияние переменной (по радиусу канала) вязкости, связанной с течением разнородных сред, на поле скоростей при ламинарных и турбулентных течениях. Рассматривались слоистые течения; было определено понятие постоянной «эффективной» вязкости, которая заменяет переменную - в определенном расчетном смысле.
Предлагаемый подход позволяет обходиться вместо достаточно сложного расчета турбулентного течения более простым расчетом ламинарного течения с флуктуациями, моделирующими вихревую засыпку. Используются уравнения Навье-Стокса без осреднения. Описанные аналитическими выражениями включения автоматически распространяются на всю область при любой ее форме. Данным методом могут рассчитываться как жидкие и газообразные (двух и многофазные) среды, так и неоднородные твердые тела (композитные, пористые и пр.). Можно учитывать изменения неоднородностей (их формы и размеров) в ходе физического процесса, т.е. ставить и решать широкий круг задач.
Для формирования расчетного поля используется понятие о характеристической функции 9(М), определенной на множестве ^ так: 9(М)=1, если точка М принадлежит множеству и 9(М)=0 в противном случае. Характеристическая функция может быть любой размерности. На расчетной сетке она образует пространство нулей с вкрапленными «островами» единиц, которые являются «носителями» определенных свойств.
Характеристические функции работают следующим образом. Зададим в качестве некоторой характеристики поля, например, вязкости д формулу
где 9 - Ц = Ц + (Ц2 -Ц,)0 ,
характеристическая функция системы включений. Из нее следует, что характеристикой основного поля, где 9=0, будет величина а характеристикой всех включений, где 9=1, будет величина ц2. Если нужно задать свойство д не только в объёмах флуктуаций, но и на их поверхностях (оболочках), эта величина представляется в виде
+ (М"2 -Ич)01 + ^302 . (10)
Здесь Д1 - характеристика внешней среды, д2 - характеристика внутренних объёмов флуктуаций, д3 - характеристика оболочек флуктуаций, 91 и 92 - характеристические функции объёмов и оболочек.
Чтобы задать в плоской области (г,г) (осевом сечении трубы) систему круговых включений с радиусом Я, расположенных в прямоугольном порядке, с расстоянием а между их центрами по оси г и Ъ по оси г, зададим характеристическую функцию этих включений в виде:
I(г, 2) = [(г - г0)2 + (2 - г0)2^ < Я(г, 2)2 .
Эта логическая функция равна единице внутри кругов и нулю вне их. Здесь
г0 = а ([г / а]+ 0,5), ^ = Ь ([г / Ь] + 0,5)
- множество координат центров кругов. Прямые скобки в данном выражении обозначают функцию «ближайшее целое слева».
Более близким к естественному (хаотическому) расположению вихрей является шахматное распределение. Оно создает большую загруженность области, чем прямоугольное. Для создания дополнительных рядов включений используем ту же функцию /, но сдвинутую относительно первой на величину т по оси г и п по оси г, а именно функцию ф(г,г) = /(г - т,2 - п). Сумма функций Дг,г) и ф(г,г) создаст шахматное распределение круговых включений. Вязкость всей структуры задается выражением
Ц = + - Ц0 )(I(Г, 2) + ф(г, 2)) .
В результате во всех точках вне круговых включений, где У(г,г) = ф(г,г)=0, получаем ц = ц0, т.е. вязкость равна вязкости текущей жидкости. В самих включениях имеют место условия У(г,г) = 1 или ф(г,г) = 1 (поскольку эти круговые включения не пересекаются), и вязкость оказывается равной заданной вязкости внутренней части круговых включений, т.е. ц=ц2. Надо отметить, что сформированные таким образом флуктуации вязкости представляют собой как бы комки вещества и двигаются вместе с потоком, но с несколько иной скоростью, получаемой в результате решения задачи.
Рассмотрим возможность вводить зависимость от координат непосредственно в характеристическую функцию. Для этого требуется изменить свойства характеристической функции, а именно, положить, что на множестве ^ 9(М) ^0 (вместо 9 (М) = 1). Примером может служить функция, представленная на рис.3, которая дает шахматное распределение синусоидальных по амплитуде пятен вязкости (или плотности) в плоской области. Здесь каждый носитель - это прямоугольник со сторонами а и Ъ. На рисунке представлены изоконтуры этого распределения при взгляде сверху. Существенно увеличив вязкость, можно придать флуктуациям свойства твердого тела.
Рассмотрим вопрос о введении в рассмотрение (при необходимости) границ флуктуаций, их оболочек. Эти границы легко описывать 5-функциями (или их аппроксимациями). Так, 5-функцию на окружности радиуса Я можно задать выражением (а/п)/(а2(Я-г)2+1) (при достаточно большом а). Можно использовать и логическую форму. Если эту же окружность считать кольцом шириной е, то ее характеристической функцией будет логическое выражение (г <Я)(г>(Я-е)).
Рис.3. Контуры равных высот флуктуаций синусоидальной формы при их шахматном расположении
(вид сверху)
В качестве примера использования вышеизложенного метода рассмотрим моделирование потока жидкости в трубе с мелкими пузырьками газа. Рассматривается ламинарное восходящее движение воды в вертикальной трубе. Расчетная область - цилиндр диаметром 0,01 и высотой 0,02 м. Шарики пара радиусом 0,5 мм в количестве 150 равномерно распределены по объёму, создавая расходное объёмное паросодержание Р=0,1. На рис.4 видна в четырех сечениях картина вязкости системы (пузырьки пара - белого цвета). Вязкость и плотность системы соответствуют формуле (10), в которой были использованы описанные выше характеристические функции в их трехмерной логической форме.
Профиль скоростей на входе - параболический. Расчеты велись в системе Сот^'о1 ЫиШркУ81с8. Результаты расчетов сравнивались с экспериментальными данными [19]. Расчеты перепадов давления и коэффициентов сопротивления потока в диапазоне чисел Рейнольдса от 100 до 1000 дали результаты, отличающиеся от экспериментальных не более чем на 10%. В результате расчета были получены также поля скоростей и несущей среды, и расположенных в нее пузырьков - включений.
Рис.4. Вязкость воды и пузырьков пара в трубе (в четырех сечениях)
Рассмотрев метод локальных флуктуаций, вернемся к модели, в которой он применяется.
Локально-сосредоточенная модель и ее реализация
Моделью с локально -сосредоточенными сопротивлениями является система флуктуаций вязкости, дающая такие же интегральные характеристики потока, что и реальная система вихрей или вихревых трубок. Иными словами, это сравнение турбулентной «вихревой засыпки» с некоторой гипотетической засыпкой, состоящей из жидких шаров другой вязкости и, в общем случае, других размеров, но дающей тот же профиль скоростей. Введем следующее определение. Назовем эквивалентной засыпку, дающую в турбулентном потоке в трубе тот же профиль скорости и то же сопротивление, что и «натуральная» вихревая засыпка, т.е. сам поток. Критерием эквивалентности будем считать равенство интегралов от сил давления, действующих на контрольную поверхность при данном значении Re. Этот признак гарантирует одинаковость динамического влияния объекта и его модели на окружающую среду.
Моделирование одиночного вихря.
Рассмотрим условия замены (в целях расчета) вращающегося подвижного вихря или вихревой трубки также подвижным, но не вращающимся объемом другой вязкости. Эквивалентность такой замены будет достигнута, в соответствии с известным методом Шварца, если вместо убранной части расчетной области будет вставлена другая часть, обладающая тем же динамическим воздействием на окружающую среду.
Процесс выбора параметров требуемой эквивалентной замены рассмотрим на примере отдельного вращающегося шарового объема жидкости, находящегося на оси трубы и движущегося вместе с жидкостью. В данном примере диаметр трубы ^=0,2 м, высота Н=0,3 м, диаметр шара ^=0,1 м, средняя скорость потока ^а=100 м/с, плотность жидкости (газа) р=1 кг/м3, динамическая вязкость ц=10-3, число Рейнольдса по диаметру шара Re=104, на входе - параболический профиль скорости. Задача решалась в системе Сот^'о1 МиШрИг^'С как осесимметричная трехмерная модель. В результате серии расчетов получено, что изменение угловой скорости вращения вихря на пять порядков (от ю=0,001 1/сек до ю=100 1/сек) мало сказывается на величине силы давления, изменяя ее не более чем на 10-12 % от среднего значения. Даже малая скорость вращения делает вращающийся
объем почти непроницаемым для внешнего потока, и степень этой непрозрачности с увеличением скорости вращения меняется мало. Вращение вихря в других направлениях почти не меняет динамическое воздействие на него. Это явление благоприятствует условиям моделирования, делая их мало зависящими от разброса скоростей вращения вихрей. Для дальнейшей работы выбирается некоторое среднее значение динамического воздействия (при данном числе Рейнольдса).
Далее возникает вопрос о том, как найти величину вязкости рассматриваемой флуктуации, которая придаст ей необходимую интегральную силу давления в данных условиях. Для определения зависимости динамического воздействия от вязкости проводится другая серия расчетов, а затем путем сравнения графически определяется искомая величина вязкости. В нашем случае она равна ц*=0,4 Па сек. Таким образом находится эквивалентная вязкость флуктуации данного радиуса при данном значении Яе для шара.
Моделирование системы вихрей
Рассмотрев получение одиночной эквивалентной флуктуации, обратимся теперь к задаче формирования системы вихрей, конкретно - применительно к круглой цилиндрической трубе [18]. Построим систему вихревых тороидальных трубок, окружающих ось трубы и расположенных в шахматном порядке. Требуется выбрать размеры (радиусы) вихревых трубок и их вязкость. Предварительно для этого надо определить, в свою очередь, а) величину и профиль эффективной («турбулентной») вязкости в трубах и б) степень влияния вязкости включений на эффективную вязкость «смеси».
Для определения параметров распределенной турбулентной вязкости применим уравнение (3). Внутреннее сопротивление в нем представлено вторым членом в скобке правой части. Это выражение является отношением у=Е/Р - градиента сил сопротивления Е к градиенту сил давления Р (см. [3]). Численное решение уравнения (3) дает возможность найти зависимость величины у от радиуса канала при различных значениях числа Re. Результаты показывают, что величина у мало отличается от единицы на большей части турбулентного ядра и резко уменьшается только при х > 0,5-0,7. Внутреннее сопротивление Е (поскольку Р постоянно) также уменьшается в переходной области, где 9-0,5, и полностью исчезает в ламинарном подслое, где 9=0. Нужно отметить, что профиль величины у достаточно хорошо соответствует известным распределениям турбулентной вязкости.
Итак, найден профиль турбулентного сопротивления Е, но не его величина. Для определения искомой величины рассмотрим вопрос о соотношении физической и турбулентной (расчетной) вязкостей, как по величине, так и по профилю. Для этого выпишем два уравнения Навье-Стокса для стационарного течения в круглой трубе. Первое уравнение - с турбулентной вязкостью ц(г) и градиентом Рь соответствующим турбулентному режиму - (ц, (г)г^')' / г = р . Второе уравнение - с физической вязкостью ц0 и градиентом Р2, соответствующим ламинарному режиму - (ц0гаТ) / г = Р . Однократно проинтегрируем (с учетом условия на оси) оба уравнения и разделим их почленно друг на
друга. Получим (в относительных величинах) соотношение Ц<(х) = 21 ^ . Используем для
У' £
/Г 1
нахождения производных у' и у/ логарифмический профиль у = — (—1пВ) в
\ 8 к
турбулентном случае и параболический у = 2(1 - х2) в ламинарном. Коэффициенты сопротивления примем такими: ^ = (1,81§Яе-1,5)-2 (П.К. Конаков) и ^=64^е. В результате получаем такую зависимость турбулентной вязкости от физической вязкости, режима течения и радиуса канала (ф(х) - парабола)
= ЛЦ0Ф( Х)= /0Жеф( Х) (11)
В соответствии с известными данными примем, что от оси и до х=0,5 отношение ц(х) /- постоянная величина, а далее оно уменьшается по параболическому закону. Это
приводит к следующему виду логической функции ф(х):
ф( х) = (х < 0,5) + 4х(1 - х)( х > 0,5) (12)
Результаты расчетов по формулам (11, 12) хорошо согласуются с известными экспериментальными данными, определяющими турбулентную вязкость. Таким образом, с помощью анализа общеизвестных уравнений для профиля скорости и коэффициента сопротивления турбулентного и ламинарного потоков получено новое простое выражение
для величины и профиля средней вязкости турбулентного потока, т.е. получен ответ на вопрос а) данного раздела.
Рассмотрим теперь вопрос б) этого раздела - вопрос о соотношении вязкости включений и средней вязкости смеси, т.е. турбулентной вязкости потока. Из-за вращения вихрей их эквивалентная вязкость значительно больше физической, что существенно увеличивает среднюю вязкость потока. В то же время вблизи стенок происходит уменьшение сопротивления, т.е. средней вязкости, за счет уменьшения размера вихрей-зерен, разрежения турбулентной «засыпки» (увеличения ее пористости). Для оценки влияния размера включений на вязкость потока была проделана серия расчетов с переменными (от расчета к расчету) размерами зерен. Сначала рассматривался предельный случай - плотная упаковка зерен некоторого радиуса R0, имеющих физическую вязкость ц2 (при вязкости жидкости Конкретная величина R0 не имеет значения, поскольку при плотной упаковке загруженность сечения не зависит от радиуса. Последовательно уменьшая радиус зерен R, получили, что отношение средней вязкости однородной смеси к вязкости зерен пропорционально ^Ж0)3. Таким образом, чтобы средняя вязкость менялась как ф(х), нужно, чтобы радиус включений уменьшался по закону
1
Я(х)= Яф( х)3 (13)
где R - принятый радиус включений вблизи оси потока.
Следующая серия расчетов была проведена с использованием зависимости (13) и переменным числом Рейнольдса - от 104 до 107. В результате была получена следующая рекомендация для вязкости зерен-флуктуаций ц2:
»0,02(^е)5'2 Л^0 (14)
Стоит отметить, что при увеличении числа Re в 10 раз расчетная вязкость зерен засыпки увеличивается в 20 раз.
Рассмотрим результаты расчетов профилей скоростей течения с заданной таким образом «засыпкой» при различных средних скоростях течения. Рассчитывался участок трубы с переходом от ламинарного режима к турбулентному. На входе в трубу задавался ламинарный профиль скорости со средним значением, соответствующим необходимому числу Re. Расчет проводился в системе Comsol Multiphуsics по схеме расчета ламинарного течения. Расстояние между центрами зерен шахматной «засыпки» составляло величину а=0,012 м. Размер зерен изменялся по радиусу канала в соответствии с выражением (13) при R=5,6•10-3 м. Вязкость «зерен» при Re=104, 105, 106, 107 была задана как 0,8, 20, 350, 8000 Пас соответственно. Были проведены расчеты течения воды при Re=104, 105, 106, 107 в вертикальной круглой трубе диаметром 0,08 м и высотой 0,2 м. Поле вязкости потока задавалось по формулам (11-14).
На рис.5 представлена средняя скорость турбулентного течения в поперечном сечении на расстоянии 0,18 м от входного сечения при вышеназванных значениях Re. Имеет место хорошее соответствие с экспериментальными данными [14]. Некоторая негладкость кривых, их квазипериодические возмущения связаны с существенной неоднородностью расчетного поля, с присутствием в расчетном поле зерен-флуктуаций вязкости, и отражают реальную суть процесса.
Таким образом, сформированная вышеописанным образом система тороидальных флуктуаций вязкости в круглой трубе, имитирующая реальную систему вихрей (вихревых колец), дает то же расчетное поле скоростей, что и в натурном эксперименте, причем во всем диапазоне турбулентных чисел Рейнольдса.
Y 1.4
1.2
05
0J5
0.4
02
Л
0 OJOQ5 OBI ОХ) 15 0В2 0J025 0J03 0Ш5 ОШ
Рис.5. Относительная скоростьy=w/wa турбулентного течения в поперечном сечении трубы при наличии «вихревой засыпки» в зависимости от радиуса r, м. Кривые соответствуют Re=104, 105, 106,
107 - сверху вниз (в левой части рисунка)
Обсуждение
В работе представлены только основные аспекты методов реализации концепции турбулентной вихревой засыпки. Круг задач, которые можно рассмотреть предлагаемыми методами, может быть расширен. Можно ввести в рассмотрение целый ряд дополнительных параметров. Параметры структуры засыпки и параметры ее характеристических функций легко сделать зависящими от времени, координат или скоростей. С помощью этих параметров засыпку можно сделать распадающейся и возникающей, колеблющейся, эллипсовидной и т.п. Могут рассчитываться неодномерные и многофазные потоки. Кроме того, в дальнейшем возможен подход, еще более близкий к естественным условиям, а именно - переход от вязких флуктуаций к вращающимся. Дальнейшее накопление данных о размерах и взаиморасположении вихрей около стенок и других препятствий позволит в будущем отказаться от опоры на профиль турбулентной вязкости.
Компьютерная реализация предлагаемых методов достаточно проста. Разница, по сравнению с другими методиками расчетов, состоит только в аналитическом задании вязкости расчетного поля. При этом отпадает необходимость задания вихрей как геометрических объектов, а также и условий на их границах. Предлагаемая методика позволяет производить расчеты турбулентных течений по программам для ламинарных, что значительно ускоряет работу, при этом действия по поддержанию устойчивости счета сводятся к минимуму, а возможности вариантных расчетов значительно увеличиваются.
Заключение
Суммируя вышеизложенное, можно утверждать, что концепция турбулентной «вихревой засыпки» позволяет:
1) получить уравнение турбулентного движения в трубах, в явном виде учитывающее распределенное внутреннее турбулентное сопротивление потока.
2) получить приближенное аналитическое решение этого уравнения, подтверждающее параболичность профиля турбулентного ядра потока в трубах и удобное в дальнейших приложениях.
3) получить численные решения, как одномерного уравнения, так и двух- и трехмерного с аналогичным внутренним сопротивлением.
4) разработать схему строения турбулентной» вихревой засыпки».
5) рассчитывать турбулентное течение на основе этой схемы, т.е. проводить численные эксперименты над физической системой, достаточно близкой к реальности и легко изменяемой.
6) выявить и использовать глубокую внутреннюю связь между турбулентными течениями и течениями в зернистом слое.
С помощью данного подхода можно проводить численные эксперименты не только применительно к течению в трубах при другой геометрии и в других условиях. Это позволит сопоставить результаты вариантных расчетов с результатами физического эксперимента и выявлять, таким образом, внутреннюю суть турбулентных процессов.
Литература
1. Falkovich G., Sreenivasan K. Lessons from hydrodynamic turbulence // Physics Today. 2006. V. 59. iss.4. pp. 43-90.
2. Baumert H.Z. Universal equations and constants of turbulent moution // Physica Scripta. 2013. V. 155. pp. 1-12.
3. Меламед Л.Э. Уравнение турбулентного движения в трубах // Письма в Журнал технической физики. 2015. T. 41. №. 24. C. 23-28.
4. Cebeci T. Turbulence Models and Their Application. Horizons Publishing Inc. Long Beach. California. 2004. 120 p.
5. Lee, M., Moser, R. D. Direct numerical simulation of turbulent channel flow up to Re=5200 // Journal of Fluid Mechanics. 2015. V. 774. pp. 395-415.
6. Аэров М.Э., Тодес О.М. Гидравлические и тепловые основы работы аппаратов со стационарным и кипящим зернистым слоем. Л.: «Химия». 1968, 510 с.
7. Заволженский М.В., Руткевич П.Б. Развитая турбулентность в трубах. М.: Институт космических исследований Российской академии наук. 2007. № 2140. 38 с.
8. Zagarola M.V., Smits A.J. Mean-flow scaling of turbulent pipe flow // Journal of Fluid Mechanics. 1998. V.373. pp. 33-79.
9. Monkewitz P. A., Nagib H. M. Large Reynolds number asymptotics of the streamwise normal stress in ZPG turbulent boundary layers // Journal of Fluid Mechanics. 2015. V.783. pp.474-503.
10. Furuichi N., Terao Y., Wada Y. et al. Friction factor and mean velocity profile for pipe flow at high Reynolds numbers // Physics of Fluids. 2015. V.27, 9(095108). pp.1-16.
11. Yongyun Hwang. Mesolayer of attached eddies in turbulent channel flow // Physical Review Fluids. 2016. V.1. 064401. pp.1-18. doi: 10.1103/PhysRevFluids.
12. Orlu R., Fiorini T., Segalini A., et all. Reynolds stress scaling in pipe flow turbulence - first results from CICLoPE // Transactions of Royal Society. 2016. V.375 (20160187). pp. 1-6.
13. Vinuesa R., Duncan R. D. Nagib H. M. Alternative interpretation of the Superpipe data and motivation for CICLoPE: The effect of decreasing viscous length scale // European Journal of Mechanics -B/Fluids. 2016.V.58. pp. 109-116.
14. Баренблатт Г.И., Корин А.Дж., Простокишин В.М. Турбулентные течения при очень больших числах Рейнольдса: уроки новых исследований // Успехи физических наук. 2014. T. 184. № 3. C. 265-272.
15. Вигдорович И.И. Описывает ли степенная формула турбулентный профиль скорости в трубе? // Успехи физических наук. 2015. Т.185. № 2. C. 213-216.
16. Баренблатт Г.И., Корин А.Дж., Простокишин В.М. К проблеме турбулентных течений в трубах при очень больших числах Рейнольдса // Успехи физических наук. 2015. Т.185. №2. C.217-220. Доступно по: https://www.ufn.ru/ru/articles/2015/2/h/. Ссылка активна на: 15 сентября 2019 г.
17. Reichardt H. Vollstandige Darstellung der turbulenten Geschwindigkeitsverteilung in glatten Leitugen // Zeitschrift fur Angewandte Mathematik und Mechanik. 1951. V.31. N.7. pp.208-219.
18. Меламед Л.Э., Филиппов Г.А. Моделирование турбулентности как «вихревой засыпки» // Известия высших учебных заведений. Проблемы энергетики. 2017. Т.19. № 9-19. C.122-132.
19. Меламед Л.Э., Филиппов Г.А. Обобщенная формула для скорости турбулентных и ламинарных течений в трубах // Известия высших учебных заведений. Проблемы энергетики. 2018. Т.20. № 7-8. C.136-146.
20.Меламед Л.Э. Метод локальных флуктуаций и моделирование неоднородных сред // Письма в Журнал технической физики. 2016. T. 42. № 19. C.31-37.
21.Кутателадзе С.С. Теплопередача и гидродинамическое сопротивление. М.: Энергоатомиздат. 1990. 366 с.
Авторы публикации
Меламед Лев Эммануилович - канд. техн. наук, ведущий научный сотрудник, АО «Интеллект» (г. Москва). E-mail: lev.melamed@yandex.ru.
Филиппов Геннадий Алексеевич - д-р техн. наук, академик РАН по Отделению энергетики, машиностроения, механики и процессов управления, РАН (г. Москва).
References
1. Falkovich G, Sreenivasan K. Lessons from hydrodynamic turbulence. Physics. Today. 2006. pp. 43-90. doi: 10.1063 / 1.2207037.
2. Baumert HZ. Universal equations and constants of turbulent motion. Physica. Scripta. 2013;55:1-12. doi:10.1088/0031-8949/2013/T155/014001.
3. Melamed LE. An equation of Turbulent Motion in Tubes. Technical Physics Letters. 2015; 41(24):23-28.
4. Cebeci T. Turbulence Models and Their Application. Horizons Publishing Inc. Long Beach. California 2004. 120 p.
5. Lee M., Moser R. D. Direct numerical simulation of turbulent channel flow up to Re=5200. Journal of Fluid Mechanics. 2015;774:395-415.
6. Aerov ME., Todes OM. Gidravlicheskie i teplovye osnovy raboty apparatov so statsionarnym i kipyashchim zernistym sloem. L . «Chemistry». 1968.510 p.
7. Zavolzhenskij MV, Rutkevich PB. Razvitaya turbulentnost' v trubah. M.: Institut kosmicheskih issledovanij Rossijskoj akademii nauk. 2007;2140:38.
8. Zagarola MV, Smits AJ. Mean-flow scaling of turbulent pipe flow. Journal of Fluid Mechanics. 1998; 373:33-79.
9. Monkewitz PA, Nagib HM. Large Reynolds number asymptotics of the streamwise normal stress in ZPG turbulent boundary layers. Journal of Fluid Mechanics. 2015;783:474-503.
10. Furuichi N, Terao Y, Wada Y. et al. Friction factor and mean velocity profile for pipe flow at high Reynolds numbers. Physics of Fluids 2015;27(9)1-16. doi: 10.1063 / 1.4930987
11. Yongyun Hwang. Mesolayer of attached eddies in turbulent channel flow. Physical Review Fluids. 2016; 1(064401):1-18. doi: 10.1103/PhysRevFluids.1.064401
12. Orlu R, Fiorini T, Segalini A. et al. Reynolds stress scaling in pipe flow turbulence-first results from CICLoPE. Transactions of Royal Society. 2016;375(20160187):1-6.
13. Vinuesa R, Duncan RD, Nagib HM. Alternative interpretation of the Superpipe data and motivation for CICLoPE: The effect of decreasing viscous length scale. European Journal of Mechanics-B/Fluids 2016;58:109-116.
14. Barenblatt GI, Korin A.Dzh, Prostokishin V.M. Turbulentnye techeniya pri ochen' bol'shih chislah Rejnol'dsa: uroki novyh issledovanij. Uspekhi fizicheskih nauk. 2014;184(3):265-272. https://ufn.ru/ru/articles/2014/3/e/
15. Vigdorovich II. Opisyvaet li stepennaya formula turbulentnyj profil' skorosti v trube? Uspekhi fizicheskih nauk. 2015;185(2):213- 216.
16. Barenblatt GI, Korin ADzh, Prostokishin VM. K probleme turbulentnyh techenij v trubah pri ochen' bol'shih chislah Rejnol'dsa. Uspekhi fizicheskih nauk. 2015;185(2):217-220. Available at: https://www.ufn.ru/ru/articles/2015/2/h/. Accessed to: 15 Aug 2019.
17. Reichardt H. Vollstandige Darstellung der turbulenten Geschwindigkeitsverteilung in glatten Leitugen. Zeitschrift fur Angewandte Mathematik und Mechanik. 1951;31(7):208-219.
18. Melamed LE, Filippov GA. Simulation of turbulence as a vortex bacfiill. Proceedings of the higher educational institutions. Energy sector problems. 2017;19(9-10):122-132.
19. Melamed LE, Filippov GA. A summarized formula for velocity of turbulent and laminar flows in pipes. Proceedings of the higher educational institutions. Energy sector problems. 2018;20(7-8):136-146.
20. Melamed LE. Method of Local Fluctuations and Simulation of Heterogeneous Media. Technical Physics Letters. 2016;42(10):990-992.
21. Kutateladze SS. Teploperedacha i gidrodinamicheskoe soprotivlenie. M: Energoatomizdat. 1990.
366 p.
Authors of the publication
Lev E. Melamed - Joint-stock company "Intelligence", Moscow, Russia. E-mail: lev.melamed@yandex.ru.
Gennady A. Filippov - Department of power, mechanical engineering, mechanics and control processes of Russian Academy of Sciences, Moscow, Russia.
Поступила в редакцию
09 сентября 2019г.
