Расчет течения вязкой жидкости в канале с учетом изменениия фазового состояния Текст научной статьи по специальности «Физика»

УДК 532.5

Вестник СПбГУ. Сер. 1. 2012. Вып. 4

РАСЧЕТ ТЕЧЕНИЯ ВЯЗКОЙ ЖИДКОСТИ В КАНАЛЕ С УЧЕТОМ ИЗМЕНЕНИИЯ ФАЗОВОГО СОСТОЯНИЯ*

А. В. Зайцев1, Е. В. Логвиненко2

1. Институт холода и биотехнологий СПбНИУ ИТМО, канд. техн. наук, доцент, zai_@inbox.ru

2. Институт холода и биотехнологий СПбНИУ ИТМО, аспирант, katena-log@mail.ru

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

В теоретическом плане — это стандартная задача механики жидкости и газа, в общем виде сформулированная, например, в [1, 2]. Обычно решение такой задачи производят после введения различных упрощений и допущений, связанных с конкретным частным рассматриваемым случаем течения. Тем не менее, в большинстве случаев аналитическое решение получить невозможно, и оно ищется с применением приближенных методов вычислительной гидродинамики, в том числе, конечно-разностных [3, 4]. Современный уровень развития вычислительных мощностей позволяет снять главное ограничение и довести решение поставленной задачи до количественных результатов. В работе в качестве численного метода используется один из наиболее гибких и универсальных методов решения сложных задач механики и теплообмена сплошных сред — конечно-разностная аппроксимация системы дифференциальных уравнений.

При компьютерном моделировании описанного процесса течения следует проанализировать большое разнообразие разностных схем. Вследствие наличия существенных нелинейностей при численном моделировании возможны серьезные проблемы сходимости и устойчивости. В результате изучения публикаций и анализа теоретических и численных исследований была предложена методика пошагового усложнения, заключающаяся в начальном решении исходной простой задачи, в простановке которой учитываются стандартные упрощения, и в дальнейшем проведении на каждом шаге повторного анализа выбранной разностной схемы, устойчивости и сходимости. Таким образом, был проделан путь от уравнения течения вязкой жидкости в трубе, инвариантного по отношению к переносам в направлении движения, и уравнения Бюргерса до параболизованного уравнения Навье—Стокса [3].

Математическая формулировка рассматриваемой задачи представляет собой следующую систему уравнений: уравнение Навье—Стокса при отсутствии касательных напряжений в текущей среде (1), уравнение неразрывности (2), уравнение энергии (3), один из возможных вариантов уравнения состояния (4). Искомыми неизвестными функциями в этой системе из четырех уравнений являются распределение продоль-

*Доклад на Международной конференции по механике «Шестые Поляховские чтения» 31 января— 3 февраля 2012 г., Санкт-Петербург, Россия.

© А. В. Зайцев, Е. В. Логвиненко, 2012

ной проекции скорости, давления, температуры и плотности по каналу — ш(х, у, х,Ь), р(х,у,х,Ь), Т(х, у, х,Ь), р(х,у, г,Ь). При использовании стандартных обозначений система уравнений имеет вид

дш дш дт дх др д{рш)

1 др 1

р дх р

д_

дх

/х-

дш дх

+

д

дш

ду ду ) ' дг дг

+

д

дш

М"

дт срр

дх

0;

дТ дТ

дт дх

др др д дт дх дх

хдТ\ д_ дх ду

хдт

ду

д_

дх

£ = гит.

(1) (2)

(3)

(4)

В дальнейшем для определённости будем рассматривать один из конкретных вариантов физической задачи — процесс начала транспортировки сжиженного природного газа (СПГ) по трубопроводу.

Рассмотрим перекачку СПГ при температуре близкой к Т по трубопроводу большой протяжённости (Ь ^ ¿). В этом процессе даже небольшие теплопритоки к прокачиваемому по трубе криопродукту на больших расстояниях от входа могут привести к получению достаточного количества теплоты для прогрева криопродукта до температуры насыщения Тя, появлению центров парообразования и нарушению гидродинамики потока. Двухфазный переход должен быть зафиксирован в результате вычислений.

Поскольку реальный процесс прокачивания СПГ по трубопроводу происходит при умеренных скоростях, в уравнении (3) пренебрегаем диссипацией энергии. Основными принятыми допущениями, требующими в дальнейшем уточнения и учета, являются отсутствие поперечных к направлению движения потока компонент скорости и давления. Тем не менее, поскольку по предварительной оценке значения Ев достигают (1, 3 — 5) • 106, учитываем турбулентность течения, вводя в соответствие с моделью пути смешения Прандтля турбулентные вязкость и теплопроводность

р = р [0, 39 (Ет — г)]2

дш дг

Рг

срр А '

где Ет —радиус трубы; г — радиальная координата.

Начальные и граничные условия вытекают из конкретной физической постановки задачи. Первоначально труба заполнена неподвижной жидкостью (шо = 0) при давлении ро = 0,1 МПа и охлаждена до температуры То. В момент времени т = 0 в трубу подаётся жидкость при температуре Тп = То. На входе в трубу (х = 0) за некоторое время тп происходит быстрый рост давления до рп и скорости до шп в соответствии с заданным законом изменения. Для определённости примем линейный закон изменения параметров во входном сечении:

Р(т) = Ро + (Рп ~Ро)

(т) =

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

ш

п

п

п

начало парообразования в жидкости, длина является неизвестной, а на конце трубы задаётся дополнительное условие, например, давление. Задача при этом может рассматриваться как маршевая по продольной координате [3].

На стенке трубы задаётся и>|г=дт = 0. В пристеночном слое в расчётах используется молекулярная вязкость и теплопроводность в ламинарном подслое, определяемом в зависимости от величины местного числа Де, а проблема потери точности в тонком пристеночном слое с большим градиентом скорости решается введением «скользящей» величины скорости на стенке [5].

Температура жидкости у стенки трубы определяется с учётом известной из условий эксплуатации трубопровода поверхностной плотности теплового потока д, подводимого к трубе извне и передаваемого к жидкости: Л (дТ/дг)г=д = д.

При переходе от дифференциалов в уравнениях (1)—(4) к конечным разностям вводим обозначения:

Г _ си г' _ си—^^с _ си Г _ си

/ = 1 = ' /¿-1 = ; /¿+1 = ;

Г _ ги Г _ ги . Г _ ги . Г _ ги

= ; /^ + 1 = ; Л-1 = ; Л+1 = /¿¿,А+1,

где обозначает одну из сеточных функций — , , или рп;^,^ ; ¿, 3 — индексы узлов пространственной сетки в поперечном сечении канала; к — индекс узлов пространственной сетки вдоль потока; п — номер временного слоя.

После многочисленных тестовых расчетов выбрана следующая разностная аппроксимация системы уравнений:

—уравнение движения (переноса импульса):

- = + (А + А + Аз )Дг, (5)

где А = — 1/р(р — рА—1)/Дг учитывает действие сил давления; А = —— — конвективный член; Аз = 1/р( (рг-1—г-1 — + ¿+1 )/Дж2 +

(^¿-1 —¿-1 — + р^+1—^+1)/Ду2 + (ра-1 —а-1 — + ра+^а+^/Д^2) —вязкостный член, учитывающий трение;

— уравнение неразрывности:

Р = Р--д^-(6)

— уравнение энергии:

Т = Т' + (В + В + Вз)Дт, (7)

где В = 1/(рср)((р — р')/Дт — —(р — рА-1)/Дг) учитывает совершаемую давлением работу; В2 = ——(Т — ТД-^/Дг — конвективная составляющая; Вз = 1/рсрх ((^-^¿-1 — 2ЛТ + Л*+1 Т1+1)/ Дж2 + (Лд--1 Т?-1 — 2ЛТ + Л^^+О/Ду2 + (Лк-1Тк-1 — 2ЛТ + Лд+1 Тд+1)/Дг^ —член, учитывающий энергию, передаваемую теплопроводностью;

— уравнение состояния [6] (для определенности выбрано уравнение Боголюбова— Майера для сжимаемой жидкости):

Для исключения потери точности на искривлённой границе используется стандартный подход — измельчение шага расчётной сетки у границы.

Целью исследования динамики начала транспортировки по трубопроводу криогенной жидкости является определение сечения падения давления до заданного предельного значения, или условий возникновения паровой фазы в потоке — момента времени и расстояния от входа в трубу до сечения начала кипения. Приведённая модель может быть также использована для решения других задач течения в трубе, например, для исследования процесса заполнения трубы жидкостью и захолажива-ния. С этой целью следует пересмотреть начальные и граничные условия.

Полученная в результате разностной аппроксимации система уравнений (5)—(8) является существенно нелинейной. Алгоритм численного решения этой системы основан на методе итераций, во многом похож на итерационно-маршевый метод [4] и отличается высоким уровнем вложенности циклов, большим количеством узлов пространственной сетки при значительной протяженности канала. В процессе алгоритмизации и численных экспериментов выработан ряд новых элементов стратегии вычислений.

1. Градиент давления рассматривается как дополнительная неизвестная функция, и к системе (5)—(8) добавляется уравнение для предварительной оценки давления в сечении на шаг вперед:

др

Рк+1 = Рк + т^ Дг. (9)

2. Уравнения неразрывности (6) и состояния (8) используются для вычисления невязки между итерациями при расчёте распределения параметров в текущем сечении. Для этого на каждой итерации сравниваются две интегральные характеристики — массовые расходы в конечно-разностном приближении

& = ^2 Р*>э ^ ДхДу, (10)

вычисляемые с использованием плотностей, определённых по уравнениям (6) и (8).

3. На численном эксперименте доказана однозначная зависимость массового расхода в любом сечении трубы от градиента давления (рисунок), что позволило построить алгоритм решения системы, основанный на итерационном поиске градиента давления методом половинного деления с проверкой результата по невязке в определении массового расхода в текущем сечении.

4. Параллельно с применением теории пути смешения для турбулентного течения в каждом узле пристеночного слоя итеративно оценивается локальное число Яв и по

нему выбирается использование турбулентных или молекулярных значений вязкости и теплопроводности.

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

Литература

1. Лойцянский Л. Г. Механика жидкости и газа. М.: Дрофа, 2003. 840 с.

2. Волландер С. В. Лекции по гидроаэромеханике / под ред. Н. Н. Поляхова. СПб.: Изд-во С. Петерб. ун-та, 2005. 300 с.

3. Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен: в 2-х т. М.: Мир, 1990. 728 с.

4. Скурин Л. И. Итерационно-маршевый метод решения задач механики жидкости и газа // Сиб. журн. вычисл. математики. 1998. Т. 1, №2. С. 171-181.

5. Патанкар С., Сполдинг Д. Тепло- и массообмен в пограничных слоях / пер. с англ. М.: Энергия, 1971. 128 с.

6. Справочник по физико-техническим основам криогеники / М. П. Малков, И. Б. Данилов, Л.Г.Зельдович, А.Б.Фрадков. М.: Энергоатомиздат, 1985. 432 с.

Статья поступила в редакцию 26 июня 2012 г.