Спросить
Войти
УДК 532.5
Вестник СПбГУ. Сер. 1. 2012. Вып. 4
РАСЧЕТ ТЕЧЕНИЯ ВЯЗКОЙ ЖИДКОСТИ В КАНАЛЕ С УЧЕТОМ ИЗМЕНЕНИИЯ ФАЗОВОГО СОСТОЯНИЯ*
А. В. Зайцев1, Е. В. Логвиненко2
Рассматривается нестационарный трехмерный процесс течения вязкой среды в канале с учетом таких сильно усложняющих поиск решения реальных условий, как несовершенство среды, зависимость теплофизических свойств в любой точке от состояния потока, влияние внутренней вязкости потока, теплоподвод извне, турбулентный режим течения, возможность возникновения фазового перехода и др.
В теоретическом плане — это стандартная задача механики жидкости и газа, в общем виде сформулированная, например, в [1, 2]. Обычно решение такой задачи производят после введения различных упрощений и допущений, связанных с конкретным частным рассматриваемым случаем течения. Тем не менее, в большинстве случаев аналитическое решение получить невозможно, и оно ищется с применением приближенных методов вычислительной гидродинамики, в том числе, конечно-разностных [3, 4]. Современный уровень развития вычислительных мощностей позволяет снять главное ограничение и довести решение поставленной задачи до количественных результатов. В работе в качестве численного метода используется один из наиболее гибких и универсальных методов решения сложных задач механики и теплообмена сплошных сред — конечно-разностная аппроксимация системы дифференциальных уравнений.
При компьютерном моделировании описанного процесса течения следует проанализировать большое разнообразие разностных схем. Вследствие наличия существенных нелинейностей при численном моделировании возможны серьезные проблемы сходимости и устойчивости. В результате изучения публикаций и анализа теоретических и численных исследований была предложена методика пошагового усложнения, заключающаяся в начальном решении исходной простой задачи, в простановке которой учитываются стандартные упрощения, и в дальнейшем проведении на каждом шаге повторного анализа выбранной разностной схемы, устойчивости и сходимости. Таким образом, был проделан путь от уравнения течения вязкой жидкости в трубе, инвариантного по отношению к переносам в направлении движения, и уравнения Бюргерса до параболизованного уравнения Навье—Стокса [3].
Математическая формулировка рассматриваемой задачи представляет собой следующую систему уравнений: уравнение Навье—Стокса при отсутствии касательных напряжений в текущей среде (1), уравнение неразрывности (2), уравнение энергии (3), один из возможных вариантов уравнения состояния (4). Искомыми неизвестными функциями в этой системе из четырех уравнений являются распределение продоль*Доклад на Международной конференции по механике «Шестые Поляховские чтения» 31 января— 3 февраля 2012 г., Санкт-Петербург, Россия.
© А. В. Зайцев, Е. В. Логвиненко, 2012
ной проекции скорости, давления, температуры и плотности по каналу — ш(х, у, х,Ь), р(х,у,х,Ь), Т(х, у, х,Ь), р(х,у, г,Ь). При использовании стандартных обозначений система уравнений имеет вид
дш дш дт дх др д{рш)
р дх р
/хдш дх
ду ду ) & дг дг
дт срр
дТ дТ
дт дх
др др д дт дх дх
хдТ\\ д_ дх ду
£ = гит.
В дальнейшем для определённости будем рассматривать один из конкретных вариантов физической задачи — процесс начала транспортировки сжиженного природного газа (СПГ) по трубопроводу.
Рассмотрим перекачку СПГ при температуре близкой к Т по трубопроводу большой протяжённости (Ь ^ ¿). В этом процессе даже небольшие теплопритоки к прокачиваемому по трубе криопродукту на больших расстояниях от входа могут привести к получению достаточного количества теплоты для прогрева криопродукта до температуры насыщения Тя, появлению центров парообразования и нарушению гидродинамики потока. Двухфазный переход должен быть зафиксирован в результате вычислений.
Поскольку реальный процесс прокачивания СПГ по трубопроводу происходит при умеренных скоростях, в уравнении (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 = Рк + т^ Дг. (9)
= ^2 Р*>э ^ ДхДу, (10)
вычисляемые с использованием плотностей, определённых по уравнениям (6) и (8).
нему выбирается использование турбулентных или молекулярных значений вязкости и теплопроводности.
Статья поступила в редакцию 26 июня 2012 г.
