Разработка алгоритма решения уравнений Навье-Стокса для течения криогенной жидкости в трубе Текст научной статьи по специальности «Математика»

УДК 531+68!

Разработка алгоритма решения уравнений Навье—Стокса для течения криогенной жидкости в трубе

Канд. техн. наук А. В. ЗАЙЦЕВ

Санкт-Петербургский государственный университет низкотемпературных и пищевых технологий

191002, Санкт-Петербург, ул. Ломоносова, 9

Finite-difference realization of a differential equation system for movement, continuity and energy has been developed jointly with equation of state to calculate three dimensional transient profiles of longitudinal velocity, pressure, temperature and density of a viscous medium flowing inside a tube. An algorithm and an appropriate computer program have been elaborated and numerical values for a concrete cryogenic substance — one-phase (liquid) methane, with specific initial and boundary conditions — have been obtained.

Keywords: Navier-Stocks equation, longitudinal velocity, pressure. Ключевые слова: уравнение Навье-Стокса, продольная скорость, давление, температура, плотность вязкой среды.

Краткая постановка обшей задачи

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

При такой постановке задача сводится к поиску решения приведенных ниже уравнений [1]:

— основного уравнения Навье—Стокса динамики вязкого газа в проекциях:

ди ди ди ди

=

от ох ду ог

= FX +

р [ dz дх

ди дх

+

+

ду

ди dv дх

д

+ 1Гг

ди dw дх

2д_

3 дх

ди dv dw \ дх ду dz J

Ь

(1)

dv dv dv dv дт дх ду dz

= F„

И

dp d

dz dx

+2— ( —\ + — ду \ ду J dz

du dv <9y дх

dv dw

^ dz dy

2 d_ 3 dy

du dv dw ^ dx dy dz

dw dw dw dw dr dx dy dz

p \ dz dx

d_ dy

dv

Мг + т

dz

dw dy

2 c)_

3 dz

du dv dw \ ^ dx dy dz J

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

др_ д{ри) d{pv) дт дх ду

уравнения энергии

d(pw) dz

ди dw\

^Yz + d^J

+

„d ( dw\

Ь

— 0;

'dT dT dT cPP I ——I- «---1- v-^—

dr

— (x—

dx V dx

dx

dy

■ w

dT

— (\—

ду V dy

(2)

(3)

^ д_ (хдТ\ +

dz \ dzJ

+qv + 2

ди\ idv\ (dw дх J \ду) \dz

+

+

dv du\ дх ду)

dw

+

ду dz)

dv\2 i du dw\2

dz дх

2 f du dv dw \' ~3 + ду + ~£h)

(5)

— уравнения состояния (одного из возможных вариантов)

- = ZRT.

(6)

Искомыми неизвестными функциями в этой системе из шести уравнений является распределение проекций скоростей, давления, температуры и плотности по каналу — и(х,у,г,т), v(x,y,z,r), w(x,y,z,r), р{х,у,г,т), T(x,y,z,r), р(х,у,г,т).

Методика решения сложных технических задач численными методами

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

Математическая постановка любой задачи в самом общем виде

Фг (/;,**) = о,

(7)

При этом констатация факта адекватности задач (7) и (8) зачастую остается без строгого научного доказательства, а решение задачи (8) необоснованно распространяется на исходную поставленную задачу (7).

Результатом решения упрощенной задачи (8) является набор значений функций }'Ахк) (поле функции), превращающих уравнения системы (8) в тождества

(9)

где Ф — функционал, определенный на множестве функций /,;

/ — неизвестная искомая функция / = /(з^); Хк — независимые переменные, которыми во многих прикладных задачах чаше всего являются время и координаты (к = 1,2,3,4);

г = 1,2,...,т — номер уравнения в общей системе уравнений, описывающих изучаемое явление; ] — 1,2,..., п — номер неизвестной функции.

Для выполнения условия однозначности решения необходимо выполнение равенства тп~ п.

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

Данный факт также следует подтвердить соответствующими вычислениями.

Когда систему не удается решить аналитически, находят численное приближенное решение с заданной точностью.

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

Д = Ф4(/;,а*). (Ю)

Если невязка оказывается в каждой точке Хк области определения достаточно малой величиной с заданной точностью, то найденный набор значений и /., является решением исходной задачи (7). В противном случае следует признать совокупность сделанных допущений ошибочной.

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

Течение вязкой жидкости в трубе, инвариантное по отношению к переносам в направлении движения

Допущения: теплофизические свойства постоянны (адиабатный процесс); поперечные (пульсационные) составляющие вектора скорости иии пренебрежимо малы по сравнению с продольной составляющей ю.

Следует определить продольные скорости в сечении, перпендикулярном направлению течения (оси г), т. е. решить уравнение

d w ц f d2w dr p \dx2

d2w d y2

1 dp

p dz

= 0.

(H)

&i(fj,xk) = 0.

(8)

Если dp/dz = 0, то в конечных разностях получаем

w.

п+1

i J

W¿

w

к и - m,j

■ w.

2w?

Ax2

ui,j+1

+

Ay 2

Ar.

(12)

Нижние индексы i, ] соответствуют узлам пространственной сетки Хг и у у, верхний индекс п — текущему узлу временной сетки. Производная по времени аппроксимирована разностью вперед (маршевая координата), вторые производные по координате — центральными разностями.

Граничные условия: слои вязкой жидкости, прилегающие к поверхности стенки канала, имеют скорость сгия, где 0 ^ с ^ 1 — коэффициент проскальзывания; иия — скорость в ядре потока. Частный случай при с = 0 опи-сываетусловие «прилипания», при с = 1 — условие «проскальзывания».

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

Тестовые расчеты произведены для течения жидкого метана.

Численный эксперимент, как и ожидалось, показал существенное влияние параметров разностной схемы на устойчивость и сходимость итерационного процесса. Так, при Ах = Ау = 0,1 м шаг Дт должен быть не больше 0,002 с. Уменьшение вязкости с 1 до 0,001 приводит к устойчивому решению при шаге по координатам 0,02 м и по времени 0,01 с.

Исследование сходимости и устойчивости

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

Решение системы разностных уравнений типа (12) производим итерационным методом, что позволит впоследствии перейти к более сложным уравнениям с различными видами нелинейностей. Считаем, что итерационный процесс сошелся, если во всех узлах разностной сетки отличие значений искомых функций на итерациях с номерами к + 1 и к не превосходит некоторой заранее заданной малой величины е, т. е. если - и? -I < е

I 3 \

для всех г,

Уравнение (12) допускает теоретическое определение условия устойчивости по методу Неймана (разложение решения в ряд Фурье) [2, 3]. Получаем

рАт

1

1

(Axf + (Ayf

1

(13)

Для достижения сходимости в сложных случаях возможно применение различных аппроксимационных формул для производных, изменение вычисляемых значений искомых функций на текущей итерации (аналог метода Зейделя) и другие численные эксперименты.

Анализ сходимости разностной схемы (12) при учете dp/dz = const проиллюстрирован на графике рис. 1.

Таким образом, численный эксперимент показывает следующее:

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

— константа в уравнении (градиент давления) не влияет на скорость сходимости;

— чем меньше диаметр трубы, тем сложнее подобрать соотношение шагов. Пусть труба имеет радиус 5 см, тогда, согласно графику, шаг по координате Ах = 1 см соответствует шагу по времени Дт = 10,5 с, но это грубая сетка по координате. При задании Ax = 1 мм получаем Дт = 0,105 с (628 итераций), т. е. расчет медленно изменяющихся нестационарных процессов может потребовать много вычислительных ресурсов.

График на рис. 2 отражает область устойчивого решения тестовой задачи.

Использование уравнения неразрывности

На следующем этапе вводим независимую искомую функцию — поле давлений. Тогда dp/dz ф const, и для нахождения двух неизвестных р и w следует рассматривать уравнение (11) совместно с еще одним уравнением, в частности с уравнением неразрывности.

Уравнение неразрывности (4) получено из закона сохранения вещества, которое для произвольного сечения для рассматриваемой задачи может быть записано в интегральном виде:

pwd F = const = Go,

(14)

где F — площадь поперечного сечения;

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

pijWijAxAy = G0

(15)

и для каждого сечения (шага по г) методом последовательных приближений находим dp/dz, превращающее уравнение (15) в тождество.

На рис. 3 приведен вариант визуализации результатов расчета поля продольных скоростей. Динамика изменения скорости в различных точках сечения потока показана на рис. 4.

Ц 700

а боо

а>

£ 500 о

5 400 т 300 200 100 о

-100

Ах = / 0,01 м

/

0,02

0,0 5

/

) 0,1

50 100 150 200 250 300 350

Шаг по времени, с

20,12

<и

х0,10

5

Ч

§0,08 у

с 0,06

и

03

Эо,04 0,02 0

,'11 111

' / / /

0 200 400 600 800 1000 1200

Шаг по времени, с

Рис. 1. Необходимое количество итераций в зависимости от соотношения шагов по координате и времени

Рис. 2. Область устойчивого решения

ООООООООО С ОООООООО^ 0 0 0 0 0 0^ ■ 0 О О OJ

¡0 0 0^

■Р

¡X. ;

оооооо©оооо<

[5 0 0 0 0 0 60 С 01 I О О О О О О О I О О О О О I 'ООО! 10001 '0 0 01 М о о I

1 в # » : а

а

оео о Ь 0 0 0 0 0 ооооо 00000000 10000 0 00000 0000000000000 0"

... * В я » * да я в я = = в да да =

я я » ш я да я « я я ® В а да в

# # К .1 # в * * в # в -4 = в в #

# # # * # # # # в ш » я я # ж

* в # да # «6 * в # « * в я = в я

# в # # # # в # # в # с я * в

« я # я » ю я в я г> » в да да в

ч в # в да <* # ш в # » # « « да да

ш « в // * » # # « # * ш * «Ш

да « * * # в # в # # Ь и я я

# в # ♦ * # в * *> # в * * :§ 4 да #

# да # в в ¡0 # # в # * # да да да ш

« * * * в в в а в # # » ш да да в

# в а * л в в # ф # # в # # в а я

# * « я а # » # # в да # # №

е » * в да # в # * * # « «г я я

в * в О * * # л в # # в я # да

* # Я Я # # *> я ё в » в в да в

* # # * в # « # # # « # я = в я в

# # да # # в * в » # в я да в да в

Щ * # « # # # * * я я в я я = да да

я я « в # » в # я я * я « в а = с

да я в я я да я я я и я в да да да

да я к я в с в к я в я да в

* да - = * * - = в вНН

8

да в я я я я я да 8

9 3 о

Л > о о <

0 О 4 ООО* 0 0 0 0 0 1 0 0 0 0 01 ООО! © 0 0 0 0 0 0! 000000000! 00000000000000!

ПУСК

Скорость

Цвет

■II

м/с 36,2 30,4 24,7 18,9 13,1 7,4 1,6

СБРОС

1,68 14,95 26,79 33 35,52 36,13 35,52 33 26,79 14,95 1,68

Поле скоростей

(г = 2 см, х = 16 мин)

Рис. 3. Поле продольных скоростей: 1, 2, 3, 4 — точки сечения потока

12 Время, мин

Рис. 4. Динамика изменения скорости: 1, 2, 3, 4 — точки сечения потока (см. рис. 3)

Параболизованное уравнение Навье—Стокса

Параболизованные уравнения Навье—Стокса применяют для внутренних течений, в которых можно выде-

лить преобладающее направление. Компонента скорости в этом основном направлении должна быть положительной, т. е. обратное течение в направлении основного потока запрещено. На компоненты скорости вторичного течения ограничений нет. Как и для всех форм па-раболизованных уравнений Навье-Стокса, диффузией в продольном направлении пренебрегают [2, 4|.

В качестве дальнейшего развития алгоритма расчета рассматриваем замену уравнения (12) параболизованным уравнением Навье—Стокса при отсутствии касательных напряжений в текучей среде:

дги дъи 1 др V (д2ы д2и> д2и>4

дт дг р дг р \ дх2 ду2 дг2,

Разностный аналог последнего:

Ь)к+1 -IV 1 Рк+1 - Рк ,

-А--Ь и>-А- =---X--Ь

Дт Аг р Аг

.'1Уг-1 - 2\у + ц>.,_1 - 2т + ги^+х

1 Дх2 + др

и;*,-! - 2ги +

дг2 ;

(17)

Здесь для упрощения записи введены следующие обозначения:

>-1.

т+1 = К+1 ,о,к\ щ-х = = <¿+1,*;

= ™к+1 = ^.к+х-

Производная по времени аппроксимируется разностью назад. Разрешая уравнение (17) относительно У)к+1, получаем маршевую задачу по координате г. Аппроксимация дш/дг разностью вперед, т. е. замена первого члена

71+1

V) - V)

уравнения (17) на —--, где V)

V)-

■ позволяет раз-

Ат -----

решить уравнение относительно г«' и получить маршевую задачу по времени. Численный эксперимент подтвердил преимущество аппроксимации дш/дг разностью вперед.

и/ = ги +

1рк+х ~Рк Р

Аг

Юк+1 - но

ю----Ь

Дг

и за счет этого — к некоторому сглаживанию поля скоростей в сечении.

Усложнение расчетной схемы приводит к усложнению компьютерной программы и, как следствие, к повышенной чувствительности с точки зрения устойчивости схемы к величине скорости на входе и соотношению шагов по времени и по координате вдоль оси г. Следует также учесть изменение определенного выше условия сходимости. С увеличением скорости потока на входе надо увеличивать шаг по г, увеличивать Ар/Аг, уменьшать шаг Дт. Возможность увеличения расчетного компьютерного времени в допустимых пределах следует учитывать при выборе расчетной схемы и соотношений шагов расчетной сетки.

Для сложных нелинейных уравнений теоретическое исследование устойчивости также вызывает определенные затруднения. Поэтому более продуктивным является применение численного эксперимента.

Уравнения неразрывности (16) и (4) используются в форме уравнения (14). Получено: сходимость обеспечивается в узком диапазоне Дт и 4Дг2. Расчеты производились в диапазоне изменения Дт от 0,16 до 1600 с. Сходимость констатировалась при выполнении в каждом расчетном сечении условия (3 = бо с заданной точностью.

Уравнение энергии. Решение системы уравнений

Рассмотрим уравнение энергии (5), пренебрегая дисси-пативной функцией, а также диффузионным переносом энергии, излучением и т. п.:

+ Р I -~--Ь —-А О----Н

+

Ах2 ' Ау2

у)к-\ -2ги + и>к+1 \

Д,г2

Дт.

(18)

'дТ дТ дТ дТ\ срР I ~—Ь и~—Н V——Ь I =

дт др

дх

др дх

ду

др ду

V)

дг ]

дрл дг

Обозначим

Аг

Тогда

1 Рк+1 ~ Рк . р Аг '

Л2 = —ы

У}к+1 - ю

Аг '

/гУг-1 — 2ш + №¿+1 — 2ги + и)т+1

Лз = Ч-А^-+-Ду2--

•шк-1 -2 гу + у)к+1 \ Аг2 ) '

V)' = ъ) + {Ах + А2 + А3)Ат.

+

А (х—\ А (х—\ —

дх \ дх) ду \ ду) дг \ дг

+ Чу

С учетом принятых ранее допущений пренебрегаем и, V, градиентом давлений по х и у, зависимостью тепло-физических свойств от состояния. Пусть qv = 0. Тогда

дТ

дТ

др дА^Т д^Г д2Т\

дт дг) \ дх2 ду2 дг2) '

Расчеты показывают, что введение в уравнение конвективного члена ъи(ди!/дг) приводит в разностной аппроксимации к преобладанию комплекса А2 наддругими

или

дТ дТ I [др др\

= +- + 1 +

дт дг срр \ дт о г,

А (д2т

срр \ дх2

д2Т д2Т^

ду2 ' дг2 ) ' Конечно-разностный аналог этого уравнения

(19)

Т'-Т Тк+1-Т | 1

—х-= -х--1--

Дг Дг срр

Р

Ат

Р . Рк+1 + и)

Аг

Теперь алгоритм решения системы уравнений (3—6) при заданных ограничениях заключается в совместном решении разностных уравнений (15), (18), (20-22). Выражение (21) не является самостоятельным независимым уравнением, оно лишь раскрывает найденное неявным образом значение др/дг через конечно-разностное приближение.

+а

Тг_1 -2Т + Тг

¿+1

Гд-1 -2Т + Т.

3 + 1

Ах2

А У2

+

+

Тк-1 — 2Т + Тк+1 \

Аг2

I 5

(20)

где Т = Т-1] к,Т' = ..., аналогично уравнению (17); а — коэффициент температуропроводности; р — давление в текущей точке (г^,к) в настоящий момент времени;

р' = р(Т') — давление в текущей точке в момент времени п + 1.

Проанализируем изменение давления вдоль трубы. На каждом участке дг имеем найденное итеративно при решении уравнений движения и неразрывности значение др/дг. В конечных разностях падение давления на участке Аг в момент времени п + 1 равно

АР=^Аг.

дг

Сумма падений давления равна суммарному гидравлическому сопротивлению. В задаче появляется ограничение ЕДр < рвх, которое определяет ограничение длины трубы при данном трении (проскальзывании) и напоре на входе.

Трехмерное поле давлений, начиная с сечения г = Аг (слой к + 1 = 2), находим из выражения

др л

Рк+1 =Рк +

дг

(21)

Появляется возможность на каждой итерации в каждом узле по известным давлению и температуре уточнить из уравнения состояния (6) плотность р'/1/1 и другие теплофизические параметры:

„п+1 = /ИТ>^ п+1 '

(22)

Выводы

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

2. Разработаны алгоритм и реализующая его компьютерная программа и получены численные результаты для конкретного криогенного вещества — метана — в однофазном (жидком) состоянии с конкретными начальными и граничными условиями.

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

Список литературы

1. Лойцянский Л. Г. Механика жидкости и газа: Учеб. для вузов. 7-е изд., испр. — М.: Дрофа, 2003.

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

3. Самарский А. А. Теория разностных схем. — М.: Наука, 1983.

4. Роун П. Вычислительная гидродинамика. — М.: Мир, 1980.