CC BY
50
УДК 532.5+664
Численное исследование математической модели
течения в канале
Зайцев А.В. [email protected]
Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики.
Институт холода и биотехнологий
Представлена математическая модель течения реальной жидкости в канале со многими усложняющими факторами. Описаны особенности алгоритма. Приведены различные результаты численного моделирования течения в трубе.
Ключевые слова: математическая модель, течение, реальная жидкость, численный эксперимент.
Numerical Research of Mathematical Model of the Current in the Channel
Zaitsev A.V. [email protected]
St. Petersburg National Research University of Information Technologies, Mechanics and Optics.
Institute of Refrigeration and Biotechnologies
Here is given the mathematical model of the current of real fluid in the channel with many complicating factors. The algorithm features are described. Different results of numerical modeling of the current in the tube are submitted for consideration.
Keywords: mathematical model, current, real fluid, numerical experiment.
Рассмотрим течение вязкой жидкости в канале. В качестве базового решения задачи примем вариант, описанный в [1]. Основными целями создания математической модели течения являются:
- исследование влияния различных усложняющих расчёт факторов: нестационарность трёхмерного процесса течения вязкой среды в канале, наличие объёмных сил и внутренних объёмных тепловыделений, несовершенство среды, зависимость теплофизических свойств в любой точке от
параметров состояния потока, влияние внутренней вязкости потока, возможные фазовые переходы, различные гидродинамические режимы течения, турбулентность, теплоподвод извне и др.;
- численные исследования условий устойчивости и сходимости при решении сложной задачи методом конечных разностей;
- разработка методов решения подобных задач для последующего моделирования и изучения на микроуровне (уравнения в частных производных) большого класса реальных технических задач.
Обычно решение таких сложных задач механики жидкости и газа в общем виде производят после введения различных упрощений и допущений, связанных с конкретным частным рассматриваемым случаем течения. Тем не менее, в большинстве случаев аналитическое решение получить не возможно и оно ищется с применением приближённых методов вычислительной гидродинамики, в первую очередь - конечно-разностных [2, 3]. При этом разрабатываются различные уникальные модели и методы, как, например, при моделировании течения неньютоновской жидкости в щелевом канале шнековой машины [4].
На данном этапе примем некоторые допущения, предполагая возможность их исключения в дальнейшем:
- отсутствие объёмных сил;
- отсутствие внутренних тепловыделений;
- в отдельных случаях расчёт будем производить для потока жидкости вплоть до наступления фазового перехода, оставив численные исследования течения двухфазной смеси на последующие этапы проведения работы;
- диссипация энергии пренебрежимо мала;
- отсутствие поперечных к направлению движения потока компонент скорости и давления.
Тогда математическая формулировка представляет собой следующую систему уравнений [5]: параболизованное уравнение Навье-Стокса (1), уравнение неразрывности (2), уравнение энергии (3), один из возможных вариантов уравнения состояния (4):
г дм дм
--у м— =
дт дх
др д(рм)
1 др V —— + —
р дх р
( л2 -,2 Л
д м д м д м
дх ду дх2
дт
срР
дх
дТ_ ~дт
■ + м-
= 0;
дТл
дх
др др д — + м— + — дт дх дх
с
X
V дх
дТ} д +
ду
\дТл
р/р = ZRT,
V
ду
+ ■
У
_д_
дх
X
дт_
дх
(1) (2)
(3)
(4)
где х, у, х, т - координаты в прямоугольной системе и время; V, р, ср, X -теплофизические свойства потока - динамическая вязкость, плотность,
теплоёмкость, теплопроводность; w (х, у, т), р (х, у, х), Т (х, у, х),
р( х, у, 2, т) - искомые неизвестные функции системы уравнений -
распределение продольной проекции скорости, давления, температуры и плотности по каналу; 2 - коэффициент сжимаемости; Я - газовая постоянная.
В результате предварительного анализа литературных источников, собственных аналитических и численных исследований была выработана методика постепенного усложнения аппроксимирующей разностной схемы, позволяющая решить серьёзные проблемы сходимости и устойчивости вследствие наличия существенных нелинейностей при решении системы разностных уравнений.
При переходе от дифференциалов в уравнениях (1)-(4) к конечным разностям вводим обозначения:
г _ гп гг _ тп—1,г _ ГП Г _ ГП
/ _ Ji, j ,к; I ~ Ji, j ,к ; Ji—l, j,к ; -1г+1~ Ji+\, j ,к;
г _ гп г _ гп . г _ гп г _ гп
Jj—1 = Ji,j—1,к; 1 j+1 = Ji, j+1,к; 1к—1 = Ji,j,к—1; 1к+1 = Ji,j,к+1,
где Inj к обозначает одну из сеточных функций - wrn j к, р"^ к, Тп к или р;" к;
г, j - индексы узлов пространственной сетки в поперечном сечении канала; к -индекс узлов пространственной сетки вдоль потока; п - номер временного слоя. В итоге в качестве одного из вариантов разностного аналога системы (1)-(4) можно привести систему (5)-(8):
w = w' + (А1 + А2 + А3) Ах;
р = р' — ри~— рк—1П—1 Ах;
А2
Т = Т' + (В + В2 + В3 )Ат; р/р = 2ЯТ,
(5)
(6)
(7)
(8)
1 р — рк—1
где А1 =---^ - определяет перенос за счёт сил давления;
р Аг
w — wk 1 А2 = — ^ к—1
Аг
- конвективный член;
Аз =
_1 — 2ци' + +1^-+1 + Ц j-lwj—1 — + Ц j+lwj+1 +
Ах2
В =
рср V
работу;
Цк—lwk—1 — + +1Щ+1 Аг 2
1 ( р — р Р — рк—1Л
Ат
— w-
Аг
Ау 2
- вязкостный член, учитывающий трение; член, учитывающий совершаемую давлением
T — т
B2 = — ^—1 - конвективная составляющая переноса энергии;
Bз
'Х^м — 2ХГ + , — + Х^+1
+ —---т----- +
V
Ах2 Ау2
■ Х к—1Гк—1 — + Х к+{гк+1Л + Аz2
- член, учитывающим энергию, передаваемую
теплопроводностью.
Формулировка начальных (НУ) и граничных (ГУ) условий предполагает однозначное подробное описание физического состояния жидкости во всех точках рассматриваемой области в начальный момент времени, задание конструктивных параметров канала (ориентация, размеры, материал и др.), обоснованное задание законов изменения параметров на границах рассчитываемой области. В то же время разные авторы формулируют НУ-ГУ для одних и тех же задач совершенно по-разному, не нарушая при этом физической сущности и математической точности получаемого решения.
Для демонстрации результатов численного эксперимента рассмотрим процесс закачки реальной жидкости для транспортировки при температуре близкой к Ts в трубопровод достаточно большой протяжённости ^ >> д). Приведённые ниже графики иллюстрируют качественные характеристики потока, полученные расчётным путём. Количественные характеристики при этом отражают влияние теплофизических свойств потока, принятых при расчёте.
Начальные условия. Первоначально труба заполнена неподвижной жидкостью = 0) при давлении pн и температуре ^ В момент времени т = 0 в трубу подаётся жидкость при температуре ^ = ^ На входе в трубу (г = 0) за некоторое время тв происходит быстрый рост давления до pв и скорости до в соответствии с заданным законом изменения. Для определённости примем при 0<х<хв линейный закон изменения параметров во входном сечении:
P (т) = Pн + (Pв — Pн)—, ™ М = %—.
Тв Тв
Граничные условия. Помимо скорости и давления на входе в трубу для однозначной разрешимости системы следует задать скорость и температуру на стенке трубы. В случае, когда необходимо определить длину трубы, на которой возникают ограничения, накладываемые на физический процесс, например давление падает ниже допустимого или в жидкости начинается парообразование, длина является неизвестной, а на конце трубы задаётся дополнительное условие, например, давление. Задача при этом может рассматриваться как маршевая по продольной координате [2].
Один из распространённых вариантов задания скорости на стенке трубы -определение её в зависимости от скорости в ядре потока w\ = =0, где
0 < с < 1 - коэффициент проскальзывания, аналогичный понятию «скользящих» величин в [6].
1
Температура жидкости у стенки трубы определяется с учётом известной поверхностной плотности теплового потока q, подводимого к трубе извне: Т|г _к _ T + q/а, где T - средняя по сечению температура потока; а -коэффициент теплоотдачи от стенки к потоку, вычисляемый по эмпирическим выражениям, типа а _ A Reй Ргс (Рг/Ргст (г/ dэкв )е.
Алгоритм численного решения сформулированной задачи является неким аналогом метода Зейделя и заключается в итеративном решении нелинейных уравнений с уточнением отдельных членов уравнений на различных этапах вычисления. Алгоритм отличается высоким уровнем вложенности циклов, большим количеством узлов пространственной сетки при значительной протяжённости канала. На рис. 1 в качестве примера приведён фрагмент блок-схемы одного из вариантов разработанной программы.
Рис. 1. Укрупнённый алгоритм ядра подпрограммы расчёта параметров потока (отражены заголовки и окончание циклов без включения тела цикла)
Отметим некоторые важные элементы алгоритма, выработанные в процессе множества численных экспериментов:
1. На численном эксперименте доказана однозначная зависимость массового расхода в любом сечении канала от градиента давления (рис. 2), что позволило построить алгоритм решения системы, основанный на итерационном поиске градиента давления методом половинного деления с проверкой результата по невязке в определении массового расхода в текущем сечении.
2. Градиент давления ф/дz рассматривается как дополнительная неизвестная функция, и к системе (5)-(8) добавляется уравнение для предварительной оценки давления в сечении на шаг вперёд:
Ф Л
Рк+1 = Рк.
С2
(9)
3. Уравнения неразрывности (6) и состояния (8) используются для вычисления невязки между итерациями при расчёте распределения параметров в текущем сечении. Для этого на каждой итерации сравниваются две интегральные характеристики - массовые расходы в конечно-разностном приближении
а=ЕРл- ЛхЛУ
(10)
вычисляемые с использованием плотностей, определённых по уравнениям (6)
и (8).
4. Учёт трения на стенке канала производится путём введения коэффициента проскальзывания 0 < с < 1. Влияние трения на параметры потока отражают графики на рис. 3.
а
с
80
40
0 400 800 _ dp кпа
dz, м
Рис. 2. Зависимость массового расхода от градиента давления
А кПа
300
200
100
с=0,5
с=0^Х ^чС= 0 7 -1,0
с=
0 200 400 600 ^ м
Рис. 3. Влияние трения на давление
5. Параллельно с применением теории пути смешения для турбулентного течения в каждом узле потока итеративно оценивается локальное число Re и по нему выбирается использование турбулентных или молекулярных значений вязкости и теплопроводности.
6. Для исключения потери устойчивости итерационного процесса вычисления при возникновении в отдельных узлах фазовых переходов со скачкообразным изменением теплофизических свойств вводится сглаживание этих свойств, например, по локальной области вокруг рассматриваемого узла.
Некоторые результаты расчётов приведены на рис. 4-7.
0
4
г = 0
г = 0,5 К
г = 0,8 К
г = R
w
8
12 т, мин
Рис. 4. Динамика изменения поля продольных скоростей у входа в канал
г=0 ,
м/с 3,2
2,8
2,4
2,0
0
30
— 1 1 200 м -
у
400 м 600 м
60
90 т, с
Рис. 5. Динамика изменения продольной скорости вдоль канала
8
Из-за разницы давлений в трубе и во входном сечении в момент подачи жидкости в трубу наблюдается скачок параметров на входе (рис. 4) и его быстрое распространение по трубе с постепенным затуханием (рис. 5). Характер поведения давления аналогичен изменению скорости. После скачка скорость плавно возрастает, а давление падает пропорционально коэффициенту проскальзывания и в зависимости от количества подводимой теплоты.
Т, К 104,4
104,0
103,6 103,2
0
50 м 200 м
100 ^300 м м
V
5
10 т, мин
Рис. 6. Динамика изменения температуры потока в различных сечениях
0
200
400
600
^ з кг/м
300 200 100
2, м
Рис. 7. Изменение параметров потока при возникновении паровой фазы
0
Температура у входного сечения незначительно колеблется, но на достаточном удалении от входа изменение может быть значительным (рис. 6). Рост температуры вследствие влияния теплопритока к трубе может быть нейтрализован снижением её за счёт падения давления.
При определённом соотношении температуры и давления в случае задания начальной температуры близкой к Т наблюдается возникновение паровой фазы
в отдельных узлах расчётной сетки, происходит скачок плотности в сотни раз, который влияет на все параметры и распространяется в соседние узлы за несколько итераций. Возникновение фазового перехода при определённых условиях может быть вызвано падением давления даже при отсутствии теплопритоков к трубе. В результате не выполняется уравнение неразрывности, физическая картина потока нарушена и кардинально нарушены условия сходимости и устойчивости. В дальнейших расчётах вместо локальных плотностей применяются их осреднённые значения, как это обычно делается при расчёте двухфазных потоков. Гидродинамическая картина потока при этом оценивается по локальным значениям плотностей. Характерные изменения усреднённой плотности и соответствующая реакция скорости потока при возникновении паровой фазы приведены на рис. 7.
Таким образом, численный эксперимент позволяет определить сечение и момент возникновения паровой фазы, исследовать закономерности течения жидкости при реальных условиях.
При начальном заполнении трубы жидкостью момент времени прохождения фронта потока через текущее сечение можно определить из выражения т= тz , где z - расстояние от входа в трубу; Лz - шаг
пространственной сетки вдоль трубы; ^ - средняя по сечению скорость потока. Анализ динамики заполнения трубопровода жидкостью показал близкую к линейной зависимость изменения расстояния от входа в трубу до фронта потока от времени на достаточно протяжённом трубопроводе. На рис. 8 приведены характерные кривые изменения скорости в различных сечениях трубопровода при его заполнении.
w, м/с
5,0 4,5 4,0 3,5 3,0
0 50 100 150 200 т, с
Рис. 8. Изменение скорости потока в различных сечениях трубопровода при его заполнении
Список литературы
1. Зайцев А.В., Пеленко Ф.В. Моделирование течения вязкой жидкости в трубе // Научный журнал СПбГУНиПТ. Серия: Процессы и аппараты пищевых производств (электронный журнал) / СПбГУНиПТ. - № 1. - Март. 2012. Режим доступа к журн.: http://processes.open-mechanics.com/ свободный.
2. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: В 2-х т. Пер. с англ. - М.: Мир, 1990.
3. Скурин Л.И. Итерационно-маршевый метод решения задач механики жидкости и газа // Сиб. журн. вычисл. математики / РАН. Сиб. отд-ние. — Новосибирск, 1998. - Т. 1, № 2. - С. 171-181.
4. Белецкий Е.В., Толчинский Ю.А. Течение неньютоновской жидкости в щелевом канале шнековой машины // Обладнання та технологи харчових виробництв: темат. зб. наук. пр. / Голов. ред. О.О. Шубш; Донец. нац. ун-т економши i торгiвлi iм. М. Туган-Барановського. - 2011. - Вип. 26. - 568 с.
5. Валландер С.В. Лекции по гидроаэромеханике. Учеб. пособие - Л.: Изд-во Ленингр. ун-та, 1978. - 296 с.
6. Патанкар С., Сполдинг Д. Тепло- и массообмен в пограничных слоях. Пер. с англ. М.: Энергия, 1971. - 128 с.
