МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ СОПРЯЖЕННОЙ ЕСТЕСТВЕННОЙ КОНВЕКЦИИ В ПАРЕ И ЖИДКОСТИ ПРИ БЕЗДРЕНАЖНОМ ХРАНЕНИИ КРИОГЕННЫХ КОМПОНЕНТОВ ТОПЛИВА Текст научной статьи по специальности «Физика»

УДК 536.25

БОТ: 10.18698/2309-3684-2020-3-4767

Математическое моделирование сопряженной естественной конвекции в паре и жидкости при бездренажном хранении криогенных компонентов топлива

© А.О. Городнов АО ГНЦ «Центр Келдыша», Москва, 125438, Россия

В работе представлен метод численного моделирования тепломассообмена при бездренажном хранении криогенных компонентов топлива в баках с учетом свободно-конвективных течений в жидкой и паровой фазах. Физико-математическая модель основана на приближении малых чисел Маха для пара и приближении Буссинеска для жидкости. Предложен способ численного моделирования начальной неоднородности температуры в паре. Адекватность заложенных в метод допущений и подходов подтверждена сравнением с экспериментальными данными по хранению азота и водорода.

Ключевые слова: бездренажное хранение, криогенный компонент топлива, естественная конвекция, приближение Буссинеска, приближение малых чисел Маха, приближение гомобаричности, численное моделирование

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

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

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

В паровой подушке бака могут реализовываться значительные температурные перепады [6], что ограничивает применение приближения Буссинеска для описания тепломассообмена [7]. Однако, малость скоростей конвективного движения позволяет воспользоваться приближением гомобаричности (приближение малых чисел Маха) [8-10], в рамках которого динамические перепады давления, связанные с движением, не учитываются в уравнениях энергии и состояния. Данный подход существенно упрощает теоретический анализ и численное моделирование по сравнению с полной постановкой системы уравнений Навье-Стокса для сжимаемого газа. Применение данной модели для конвекции в газе также показало наличие значительной по величине вертикальной температурной стратификации даже в случае интенсивной ламинарной конвекции водородного пара с числами Рэлея порядка 106—108 [11], причем полученные в расчетах перепады температуры значительно превышали наблюдавшиеся в экспериментах по бездренажному хранению.

В последнее время были опубликованы работы, в которых для моделирования модельных экспериментов и натурных испытаний по исследованию бездренажного хранения использовались современные вычислительные программные комплексы, позволяющие решать задачу в сопряженной многомерной постановке с учетом реальной геометрии бака [12-14]. Однако, несмотря на хорошее совпадение с данными экспериментов по давлению, в данных работах наблюдалось значительное расхождение по температурному расслоению в заполненной газом части объема бака. Необходимо отметить, что одним из факторов, учет которого может существенно влиять на свободно-конвективный тепломассообмен в паровой подушке, является стенка бака.

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

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

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

5

к

<и 03 £ Л и о

Рис. 1. Схема рассматриваемой сопряженной задачи

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

Так как в рассматриваемой постановке скорости движения газа в паровой подушке существенно меньше скорости звука, математическая модель для пара основана на приближении гомобаричности. Важной особенностью данного подхода является разделение давления на «термодинамическую», зависящую только от времени составляющую, и «динамическую», переменную как по времени, так и по пространству:

Р = РТ ^) + р (г, I).

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

Уравнения неразрывности, движения и энергии для описания процессов в паре запишем в цилиндрических координатах для осесиммет-ричного случая [15] в безразмерном виде. Для этого введем следующие соотношения для масштабов и основных критериев подобия:

Р =

Р_ Ро

Р Р Т Т и М Кй' , л (Г, г)

Р

Т

ао/ Я

я

? =

я2/ а

А =

Р =

о

—Т

Роао я

Рг =

Що с

о" Р о

Л)

ЛТо =

Ло

ср =

О)

о

Т

п

о

ЯаУ = ——яР0; ао =

аощоТо

Ло

Рос

о Р о

1 Л- Iй

Ло що

Здесь «волна» над переменными относится к размерным величинам, индекс «0» обозначает некоторые реперные значения. С учетом введенных обозначений и, опуская «волну» над переменными, уравнения для пара можно записать в следующем виде:

о 1 5

+рг^-т-

г дг

др ^ 1 Огрй ^ дрм _ д? г дг дг

дрй 1 дгрйй дрйм _ др

д? г дг дг дг

( 4 дй 2 й 2 дм I 3 дг 3 г 3 дг ,

+ Рг -

дг

(дм дй

щ —+—

\дг дг,

(2)

(3)

дрм ^ 1 дгрйм ^ дрмм _ Ф+рг д

д? г дг

дг

дг

дг

(4 2 дй_ 2 иЛ

I 3 дг 3 дг 3 г ,

1д

+

г дг

( дй дм4 ^ дг дг ,

+ (1 _р) Яа^_

(4)

Будем рассматривать пар как совершенный газ, учитывая температурные зависимости теплоемкости, вязкости и теплопроводности:

Рт =рТ.

(5)

Уравнение энергии пара запишем в следующем виде:

дрТ дг

1 дгрыТ дрмТ

дг

дг

дz

1 д

г дг

дг

1

дт

йРт йг

(6)

Уравнения неразрывности, движения и переноса тепла в жидкости запишем в тех же безразмерных масштабах:

1 дги дм Л

--+ —-0;

г дг дг

ди 1 дгии дим т.р др Кл

— +--+-= -КрЬ — + РгЬ V-

дг г дг дг УЬ дг ЬКЕ

л

д

г

дг

дм 1 дгим дмм др „

- +---— + ^— = -Крр-^ + Рг,

пдР , КУЬ

дг г дг

дг

дг

ЯаГ Рг К

К

УЬ

1 дги Л + д 2и г дг ) дг2 2,

1 д ( дм Л д2 м

--1 г— 1 + —-

г дг \ дг ) дг

л 3 УЬ

Е 2 УЬ

КЕ

дТ дг

■ +

1 дгиТ дмТ --+-

г дг дг

Кь

— т^л =

(т -1);

С1 ы.{гСт. лл

г дг ^ дг ) дг ^ дг _

(7)

(8)

(9) (10)

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

Рг, ; ^ = ВРьЪЛ рь СЬ .

1

ЛьМьТо

Ь ■ ЪТР — р0 ■

КУЬ = ;

рЬ

1

К л — ' Ь ■ 1гЕ _

КУЬ = г, ; КУЬ =

Л

рЬСЬ

р0СР 0

Перенос тепла в стенке бака будем моделировать уравнением теплопроводности [16]:

КЕ п дТ- К л

дг

1 Ц^гТ л+а^ЫТ л

г д I дг ) дг I дг ,

(11)

Здесь:

уЕ СШ0рш . у-л Лт0 .

КУШ ; КУШ 2 ;

СР 0р0 л

г — ■ 2 — ^

"Ш 0

лш 0

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

5 = 4; Б = ^; Бу = Н-; Ф = Нь

я ь я у я ' Н + Ну

На оси симметри в паре и жидкости зададим следующие соотношения для температуры и компонент скорости:

дТ_

дг

=0; й 0 = 0; ^

г-о 1г=0 ; дг

= 0. (12)

г О

На внешней поверхности стенки в качестве граничного условия зададим раномерный нагрев:

—

ОГ

- А. КГЛ ; ^П

г=1+5 сЪ

- А. (13)

г=БГ +Б1 +5,г=0

На поверхности раздела газ-стенка будем считать температуру и тепловой поток непрерывными. Граничные условия для скорости на внутренней поверхности стенки, контактирующей с паром:

й| = 0; = 0; й| = 0; = 0. (14)

1г=1 ' I г=1 ' \г=Бу +БЬ ' \г=Бу +БЬ 4 '

Постановка граничных условий на омываемой жидкостью части стенки осложняется возможностью реализации режимов теплообмена, при которых температура на поверхности раздела жидкость-стенка вырастет выше температуры насыщения [17]. В таком случае при достаточно больших величинах перегрева на этой поверхности может начаться кипение жидкости [18]. Для учета данного эффекта предложена упрощенная модель фазовых переходов. Обозначим как Тб температуру в некоторой точке на поверхности раздела жидкость-стенка. В рамках данной модели при достижении температурой Тб уровня температуры насыщения начинается кипение. Температура Тб остается равной температуре по кривой насыщения Тз плюс некоторая надбавка А Тмз, внесение которой учитывает возможность нахождения жидкости в перегретом состоянии, а парообразование рассчитывается по разнице тепловых потоков из стенки в поверхность раздела жидкость-стенка и из поверхности раздела жидкость-стенка в жидкость. Интеграл данной массы пара автоматически переносится через площадь поверхности раздела жидкость-газ. Влияние пузырей пара на течение и теплообмен в жидкости в рамках данной модели не рассматривается. Тогда граничные условия на поверхности раздела стенка-жидкость и поверхности раздела жидкость-пар запишутся следующим образом:

T < T + AT

TD < TS + AT MS

K

VL

r or}

уд r j

= K

VW

DL

\ дТ Л

AWT"

. dr j

WD

(15)

fxdt

У

dz

- к

VL

VS

a—

dz

= KVL rVL Ps WS

SL

T > T + AT

TD > TS + ATMS

к

r dT\

VL

Уд r J

DL

f д)Т~

+ 4lA = KVw AW^

У dr.

WD

(16)

fA?t

У

dz

- к

VL

VS

AdT.

dz

+ Ql = KVL rVL PS WS

SL

Bl +S

Ql = 2a j qldz.

(17)

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

Численный алгоритм решения задачи. Основу предложенного алгоритма составляет метод SIMPLE [19] для расчета переменных скорость-давление. Применение данного алгоритма к расчету параметров пара осложняется необходимостью совместного решения уравнений для скорости и уравнения Пуассона для поправок давления вместе с уравнениями для температуры, плотности и термодинамического давления пара. Данная особенность может существенно влиять на скорость сходимости решения.

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

Yt

-P

дТ u drT

dt

r

д '

dr

-w-

dT 1 - Yo dz

P VU-

1_

r dr

dT

rA

у dr

d f aT

dz

dz

(18)

Здесь:

s

Уи = 1 ЁШ + ^ ; ут=_^0_,

г дг дг 1 + /0 (-1)

Для коррекции дивергенции скорости в правой части уравнения используется следующее соотношение:

чи =

1 д ( ,дТ\ г А-

V дг)

г дг

д + —

дг

хдТ

дг

)

ут 4

(19)

Дискретизация расчетной области проводится методом контрольных объемов на разнесенной сетке. Для аппроксимации временных производных использовалась схема Эйлера 1 -го порядка, дискретизация конвективных членов проводилась с помощью схемы 2-го порядка точности с разностью против потока Самарского [21], значения диффузионных потоков на гранях рассчитывались с использованием центральных разностей. Для расчета уравнений для температуры и компонент скорости применялся метод переменных направлений, для решения уравнений Пуассона для поправок давления применялся метод сопряженных градиентов [22].

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

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

М'

Р Ш+1 _ _М

РТ = 1 Бу + Б1 +5

1 Бу +Б1 + 5

I 1

0 Бг +5

гйгйг

ш +1

; мш =| I р'гйЫг.

(20)

0 Б, + 5

Далее рассчитывается поле плотности пара из уравнения состояния. Из кривого насыщения вычисляется новая температура поверхности раздела фаз:

Тш+1 = у (р'1).

(21)

После этого, используя соотношения на поверхностях раздела жидкость-стенка и жидкость-пар, находится скорость испарения жидкости на зеркале. После нахождения скорости, из объемного интеграла

уравнения неразрывности находится новая масса пара:

мш+1 = мп +| р wsгdгdz. (22)

0

Далее проверяется сходимость температуры в емкости. В данной работе был выбран следующий критерий сходимости:

l+S By + BL +2S

*=! i

0 0

rj-im+l rj-im

jn

^Zmax • (23)

Если соотношение (23) не выполняется, то значения переменных на предыдущей итерации обновляются рассчитанными и выполняется переход обратно к расчету температуры. При выполнении (23) происходит расчет динамических уравнений методом SIMPLE c учетом рассчитанных температур, плотностей и скорости испарения на поверхности раздела жидкость-пар:

P = const = PT (0); (24)

Ts (0) = f (Pt (0)); (25)

TL=B+s= TS • (26)

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

Валидация математической модели и численного метода. Для тестирования задачи были проведены расчеты бездренажного хранения водорода и азота, результаты которых сравнивались с данными экспериментов. Первый тестовый пример — бездренажное хранение водорода в малом модельном цилиндрическом баке радиуса 0,05 м с полусферическими днищами [24] (рис. 2).

Будем рассматривать приближение реальной геометрии бака цилиндром с такой же высотой и радиусом. Толщина стенки бака — 1 мм, материал стенки — сталь 12Х18Н10Т. Теплофизические свойства паров водорода и материала стенки в безразмерном виде задавались следующими полиномами, обобщающими данные [25-27]:

Л = l - 0,000427 (T -1)3 + 0,00091 (T -1)2 + 0,85236 (T -1);

ц = 1 - 0,0018 (Т-I)3 + 0,05376 (Т -1)2 + 0,91534 (Т -1); сР -1; Т < 2;

сР --0,0027• Т3 + 0,0510• Т2 -0,1803• Т +1,6781;

70 -1

2 < Т < 10; сР - сР (10); Т > 10; Т -1 + 0,001222(Р-1)3 -0,019799(Р-1)2 + 0,149263(Р-1); ^ -1 -1,19601 (Г -1)3 -0,31561(Т -1)2 -0,21121 (г -1); ^ -1 + 0,0194(Т-1)3 -0,11004(Т-1)2 + 0,56324(Т-1); -1 + 0,14149(Т-1)3 + 0,38151(Т-1)2 + 3,30122(Т-1). Значения определяющих параметров брались следующими:

А -12,66; А -12,66; Яау- 109; Яаь -109; Ргк- 1,0; ^ - 0,15; Рг£ -1,1; 70 -1,667; К^ - 0,0172; К^ - 8,0; ф-0,75; К^ - 53,3; К^ - 2,2; Кхш - 234; К^ - 7,08; Яау-109; Яаь - 109; Ргу-1,0; ^ - 0,15; Рг, -1,1; 70 -1,667; К^ - 0,0172; К^ - 8,0; ф-0,75; КЕ, - 53,3; К^ - 2,2; К"ш - 234; КЕШ - 7,08.

Рис. 2. Схема бака из эксперимента по хранению водорода

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

Рассматривались три сетки: грубая равномерная сетка с шагом по г и г: 0,025 на 6400 контрольных объемов в паре и жидкости, сетка со сгущением у границ от максимального шага по г и г: 0,025 до

минимального шага 0,00625 размером 13298 узлов; сетка со сгущением у границ от максимального шага по г и г: 0,01 до минимального шага 0,0025 размером 48400 узлов. Стенка для всех сеток делилась на 5 контрольных объемов по толщине. На рис. 3 представлены результаты расчета роста термодинамического давления на различных сетках. Расхождение между результатами расчета давления на всех типах сеток составило менее 1 % от величины начального давления.

Рт 4

1

0 0,05 0,10 0,15 I

Рис. 3. Расчет давления

— — подробная сетка, □ — сетка со сгущением узлов, ▲ — грубая сетка

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

2 - 0,5ВЬ + 5, г - 0,5ВК + Вь + 3.

400 300 200 100 0

-100

0 0,2 0,4 0,6 0,8 1,0 г

25 20 15 10 5 0 -5

0 0,2 0,4 0,6 0,8 1,0 г

б

Рис. 4. Профили вертикальной компоненты скорости в сечении

г - 0,5В + В +5, а — г - 0,015, б — г - 0,15:

- — подробная сетка, □ — сетка со сгущением узлов, ▲ — грубая сетка

3

2

Ц

Ц

а

500 400 300 200 100

-100

0 0,2

0,4 0,6 а

0,8

1,0 г

250

150

50

-50

0,2

0,4

0,6

0,8

1,0 г

б

Рис. 5. Профили вертикальной компоненты скорости в сечении

г = 0,5^ + 3, а — г = 0,015, б — г = 0,15 : — — подробная сетка, — — сетка со сгущением узлов, — — грубая сетка

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

Для валидации модели была проведена серия расчетов задачи при различной степени заполнения. Результаты сравнения расчетных и экспериментальных данных по давлению представлены на рис. 8. Результаты сравнения расчетных и экспериментальных данных по температуре на оси симметрии бака представлены на рис. 9. Максимальное отклонение от данных эксперимента по приросту давления составило 25 %. Максимальное отклонение по приросту температуры не превысило 10 %. Существенное отклонение по давлению объясняется упрощением реальной геометрии бака в расчетах. Однако, с помощью предложенного метода, в отличии от упрощенных методик [4, 5],

w

0

w

0

получилось корректно воспроизвести эффект понижения скорости роста давления с уменьшением степени заполнения.

Рис. 6. Изолинии температуры в паре, рассчитанные на различных сетках в момент времени t = 0,15: — — подробная сетка, — — сетка со сгущением узлов, — — грубая сетка

4 3,5 3

2,5 2 1,5 1

0 0,05 0,10 0,15 /

Рис. 7. Сравнение расчетных и экспериментальных данных по росту давления:

Ф = 0,86 : — — расчет, ■ — эксперимент; Ф = 0,62 : — — расчет, ■ — эксперимент; Ф = 0,37 : — — расчет, ■ — эксперимент

0

Г

1

1,0 1,5 2,0 2,5 3,0

Рис. 8. Сравнение расчетных и экспериментальных данных по температуре, ф = 0,80:

— — расчет;

х _

эксперимент

Второй тестовый пример — бездренажное хранение азота в стальном цилиндрическом баке [23]. Схема геометрии данного бака

4

3

2

1

представлена на рис. 9. Толщина стенки бака бралась равной 1 мм. Теплоемкость газообразного азота при низких температурах хорошо описывается моделью идеального газа, поэтому в расчетах величина ср бралась постоянной и равной 3,. Зависимость теплопроводности пара в безразмерном виде аппроксимировалась следующим полиномом:

1-1-0,02733 {Т -1)3 + 0,07773 (Г -1)2 + 0,82087 {т -1), Г = 4-, Т* - 77,35К, Л = 4", Г-0,00766-Вт

Т А м • К

Зависимости теплоемкости и теплопроводности стенки от температуры аппроксимировалась следующими соотношениями:

4=1 + 0,09423 (Т -1)3 - 0,22645 (т -\)2 + 0,49988 (Г -1), Т = у, Г = 77,35К,

Лг = "Г^' Лг = 2 -

а м • к

1 = 0,00694 (г -1)3 - 0,25380 (Г -1)2 +1,08772 (Г -1)

С]Г = 1 = 0,00694^7 -1| -0,2538017 -1) +1 т = у, Т* = 77,35К,

с с *=202-Дж

т- ; уу ТЛ

% кг-К

Значения определяющих параметров брались следующими: Pr 0,75; Pri = 2,1; у = 1,4; К^ = 0,00546; КА = 17,75;

К^ = 344,4; К^ = 2,459; К^ = 1070,5;

К^ = 346,0; Яа^ = 0,6.

В эксперименте рассматривались степени заполнения бака жидкостью 70 %, 50 % и 30 %. Интегральный тепловой поток ^ в бак измерялся по расходу пара при открытом дренажном клапане на стационарном режиме. Для максимального заполнения бака он составлял 2,5 Вт, для степени заполнения 50 % — 1,2 Вт, для минимального уровня — 1,0 Вт. Необходимый для расчетов удельный тепловой поток пересчитывался путем деления интегрального теплового потока общую площадь стенки бака. В таблице 1 приведены значения

параметра А и чисел Рэлея для пара и жидкости для рассматриваемых опытов, используемые в расчетах.

й = 0,1 м 1

г--*!

Рис. 9. Схема азотного бака

Таблица 1

Значения тепловых потоков и чисел Рэлея в расчетах

Параметр 30 % 50 % 70 %

0,3 0,5 0,7

йш, Вт 1,0 1,2 2,5

, Вт/м2К 5,13 6,16 12,84

А 0,87 1,04 2,17

Яаг 5,3-107 6,3-107 1,3108

Яа £ 1,9108 2,3-108 4,7-108

В ходе экспериментов на момент начала хранения в паре реализо-вывались существенные градиенты температуры. Поэтому расчеты проводились для двух типов начальных условий для параметров пара. В первом случае задавалась однородная температура, равная температуре насыщения при начальном давлении. Во втором случае градиент температуры, реализовывавшийся в эксперименте, моделировался численно описанным ранее способом. Начальные значения температуры поверхности раздела вычислялись из кривой насыщения азота [27]. На рис. 10 представлен пример одного из рассчитанных начальных профилей температуры.

На рис. 11 показано сравнения расчетных данных по росту давления с экспериментом.

г

Рис. 10. Сравнение расчетного начального профиля температуры пара с опытными

данными р = 0,5, Q = 1,2 Вт

— — расчет; х — эксперимент

а б

Рис. 11. Сравнение расчетных (сплошные линии) и экспериментальных (точки) данных по росту давления в азотном баке, а — неоднородные начальные условия, б — однородные:

--р= 0,3 ,--р= 0,5 ,--р= 0,7

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

Из представленных данных видно, что расчётные и экспериментальные профили температуры в баке достаточно хорошо согласуются между собой. Максимальное расхождение рассчитанных значений температуры по сравнению с экспериментальными данными составило 20%.

Рис. 12. Сравнение расчетных и экспериментальных данных по температуре в азотном баке в момент времени t = 0,6, а — (р = 0,3, б — (р = 0,5:

— — расчет; х — эксперимент

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

ЛИТЕРАТУРА

[1] Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Том IV. Гидродинамика. Москва, Наука, 1986, 736 с.

[2] Полежаев В.И., Черкасов С.Г. Нестационарная тепловая конвекция в цилиндрическом сосуде при боковом подводе тепла. Известия Академии наук СССР. Механика жидкости и газа, 1983, № 4, с. 148-157.

[3] Моиссева Л.А., Черкасов С.Г. Математическое моделирование естественной конвекции в вертикальном цилиндрическом баке при знакопеременном распределении теплового потока на стенке. Известия Российской академии наук. Механика жидкости и газа, 1996, № 2, с. 66-72.

[4] Черкасов С.Г., Миронов В.В., Миронова Н.А., Моисеева Л.А. Метод расчета скорости роста давления при бездренажном хранении жидкого водорода в емкостях. Известия Российской академии наук. Энергетика, 2012, № 4, с. 151-159.

[5] Амирханян Н.В., Черкасов С.Г. Теоретический анализ и методика расчета теплофизических процессов, протекающих в криогенной емкости в режиме бездренажного хранения. Теплофизика высоких температур, 2001, т. 39, № 6, с. 970-976.

[6] Ward W.D., et al. Evaluation of AS-203 Low-Gravity Orbital Experiment. NASA CR 94045. Chrysler Corp. Space Div. Technical Report BB-3.4.3-5-101, 1967.

[7] Gray D.G., Giorgini A. The validity of the boussinesq approximation for liquids and gases. International Journal of Heat and Mass Transfer, 1976, vol. 19, pp. 545-551.

[8] Paolucci S. Filtering of Soundfrom the Navier-Stokes Equations. Sandia National Laboratories, Livermore, 1982.

[9] Лапин Ю.В., Стрелец М.Х. Внутренние течения газовых смесей. Москва, Наука, 1989, 368 с.

[10] Черкасов С.Г. О некоторых особенностях описания тепловых и динамических процессов в газах в приближении гомобаричности. Теплофизика высоких температур, 2010, т. 48, № 3, с. 444-448.

[11] Черкасов С.Г., Лаптев И.В., Ананьев А.В., Городнов А.О. Рост давления при нестационарной естественной конвекции паров водорода в вертикальном цилиндрическом сосуде с постоянной температурой нижней границы. Тепловые процессы в технике, 2019, т. 11, № 5, с. 203-215.

[12] Kartuzova O., Kassemi M., Agui J., Moder J. Self-pressurization and spray cooling simulations of the multipurpose hydrogen test bed (MHTB) ground-based experiment. Paper presented at the 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference 2014, 2014. DOI:10.2514/6.2014-3578

[13] Panzarella C.H., Kassemi M. Self-pressurization of large cryogenic tank in space. Journal of Spacecraft and Rockets, 2005, vol. 42, no. 2, pp. 299-308.

[14] Grayson G., Lopez A., Chandler F., Hastings L., Hedayat A., Brethour J. CFD modeling of helium pressurant effects on cryogenic tank pressure rise rates in normal gravity. Collection of Technical Papers — 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference, 2007, vol. 5, pp. 4956-4965.

[15] Ferziger J.H., Peric M. Computational methods for fluid dynamics. New York, Springer, 2002, 423 p.

[16] Лыков А.В. Теория теплопроводности. Москва, Высшая школа, 1967, 600 с.

[17] Вальциферов Ю.В., Дронов В.П. Численное моделирование конвективного теплообмена в тонкостенном цилиндрическом сосуде с полусферическими днищами, полностью заполненном жидкостью. Известия Академии наук СССР. Механика жидкости и газа, 1984, № 5, с. 204-207.

[18] Веркин Б.И., Кириченко Ю.А., Русанов К.В. Теплообмен при кипении криогенных жидкостей. Киев, Наук. думка, 1987, 264 с.

[19] Patankar S. Numerical heat transfer andfluidflow. Washington, DC, Hemisphere Publishing Corporation, 1980, 210 p.

[20] Quazzani J., Garrabos Y. A new numerical algorithm for low Mach number supercritical fluid. Preprint Elsevier, 2007, 10 p.

[21] Самарский А.А. Теория разностных схем. Москва, Наука, 1989, 616 с.

[22] Chung T.J. Computational fluid dynamics. Cambridge, Cambridge University press, 2010, 1012 p.

[23] Seo M., Jeong S. Analysis of self-pressurization phenomenon of cryogenic fluid storage tank with thermal diffusion model. Cryogenics, 2010, vol. 50, no. 9, pp. 549-555.

[24] Belyayev A.Yu., Ivanov A.V., Egorov S.D., Voyteshonok V.S., Mironov V.M. Pathways to solve the problem of cryogenic rocket propellant long storage in space. Proceedings of the International Aerospace Congress, Moscow, Russia, 1994, vol. 1, pp. 558-562.

[25] Григорьев И.С., Мелихов Е.З. Физические величины. Москва, Энергоатомиз-дат, 1991, 1232 с.

[26] Солнцев Ю.П., Ермаков Б.С., Слепцов О.И. Материалы для низких и криогенных температур: энциклопедический справочник. Санкт-Петербург, Химиздат, 2008, 768 с.

[27] Варгафтик Н.Б. Справочник по теплофизическим свойствам газов и жидкостей. Москва, Наука, 1972, 720 с.

Статья поступила в редакцию 03.12.2020

Ссылку на эту статью просим оформлять следующим образом:

Городнов А.О. Математическое моделирование сопряженной естественной конвекции в паре и жидкости при бездренажном хранении криогенных компонентов топлива. Математическое моделирование и численные методы, 2020, № 3, с. 47-67.

Городнов Анатолий Олегович — научный сотрудник АО ГНЦ «Центр Келдыша», область научных интересов: механика жидкости и газа, теплофизика, численные методы, математическое моделирование, криогенная техника e-mail: an.ol.gorodnov@gmail.com

Mathematical modelling of conjugate natural convection in vapor and liquid during ventless storage of cryogenic fuel components

© A.O. Gorodnov SSC Keldysh Research Center, Moscow, 125438, Russia

Method for numerical modelling of natural convection heat and mass transfer in liquid and vapor phases during ventless storage of cryogenic fuel components in tanks is presented in this paper. The proposed mathematical model is based on low Much number approach for vapor phase and Boussinesq approach for liquid phase. Method for numerical modelling of inhomogeneous initial conditions in vapor phase is given. Presented model and method is validated using experimental data for ventless storage of liquid hydrogen and nitrogen.

Keywords: ventless storage, cryogenic fuel component, natural convection, Boussinesq approach, low Much number approach, homobarisity approach, numerical modelling

REFERENCES

[1] Landau L.D., Lifshic E.M. Teoreticheskaya fizika. Tom IV. Gidrodinamika [Theoretical Physics. Vol. IV. Hydrodynamics]. Moscow, Nauka Publ., 1986, 736 p.

[2] Polezhaev V.I., Cherkasov S.G. Nestacionarnaya teplovaya konvekciya v cilindricheskom sosude pri bokovom podvode tepla [Nonstationary thermal convection in a cylindrical vessel with a side heat supply]. Izvestiya Akademii nauk SSSR. Mekhanika zhidkosti i gaza [Proceedings of the USSR Academy of Sciences. Fluid and Gas Mechanics], 1983, no. 4, pp. 148-157.

[3] Moisseva L.A., Cherkasov S.G. Matematicheskoe modelirovanie estestven-noj kon-vekcii v vertikal'nom cilindricheskom bake pri znakoperemennom raspredelenii teplovogo potoka na stenke [Mathematical modeling of natural convection in a vertical cylindrical tank with alternating heat flow distribution on the wall]. Izvestiya Rossijskoj akademii nauk. Mekhanika zhidkosti i gaza [Proceedings of the Russian Academy of Sciences. Mechanics of Liquid and Gas], 1996, no. 2, pp. 66-72.

[4] Cherkasov S.G., Mironov V.V., Mironova N.A., Moiseeva L.A. Method of calculation of pressure velocity growth at non-drainage storage of liquid hydrogen in enclosures. Izvestiya Rossijskoj akademii nauk. Energetika [Proceedings of the Russian Academy of Sciences. Mechanics of Liquid and Gas], 2012, no. 4, pp. 151-159.

[5] Amirkhanyan N.V., Cherkasov S.G. Theoretical analysis and procedure for the calculation of thermophysical processes occurring in a cryogenic vessel under conditions of nonvented storage. High Temperature, 2001, vol. 39, no. 6, pp. 905-911.

[6] Ward W.D., et al. Evaluation of AS-203 Low-Gravity Orbital Experiment. NASA CR 94045. Chrysler Corp. Space Div. Technical Report BB-3.4.3-5-101, 1967.

[7] Gray D.G., Giorgini A. The validity of the boussinesq approximation for liquids and gases. International Journal of Heat and Mass Transfer, 1976, vol. 19, pp. 545-551.

[8] Paolucci S. Filtering of Soundfrom the Navier-Stokes Equations. Sandia National Laboratories, Livermore, 1982.

[9] Lapin Yu.V., Strelec M.H. Vnutrennie techeniya gazovyh smesej [Internal flows of gas mixtures]. Moscow, Nauka Publ., 1989, 368 p.

[10] Cherkasov S.G. Some special features of description of thermal and dynamic processes in gases in the approximation of homobaricity. High Temperature, 2010, vol. 48, no. 3, pp. 422-426.

[11] Cherkasov S.G., Laptev I.V., Ananyev A.V., Gorodnov A.O. Pressure rise during unsteady natural convection of hydrogen vapor in a vertical cylinder with isothermal bottom boundary. Teplovye processy v tekhnike [Thermal processes in engineering], 2019, vol. 11, no. 5, pp. 203-215.

[12] Kartuzova O., Kassemi M., Agui J., Moder J. Self-pressurization and spray cooling simulations of the multipurpose hydrogen test bed (MHTB) ground-based experiment. 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference 2014, 2014. DOI:10.2514/6.2014-3578

[13] Panzarella C.H., Kassemi M. Self-pressurization of large cryogenic tank in space. Journal of Spacecraft and Rockets, 2005, vol. 42, no. 2, pp. 299-308.

[14] Grayson G., Lopez A., Chandler F., Hastings L., Hedayat A., Brethour J. CFD modeling of helium pressurant effects on cryogenic tank pressure rise rates in normal gravity. Collection of Technical Papers — 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference, 2007, vol. 5, pp. 4956-4965.

[15] Ferziger J.H., Peric M. Computational methods for fluid dynamics. New York, Springer, 2002, 423 p.

[16] Lykov A.V. Teoriya teploprovodnosti [Theory of thermal conductivity]. Moscow, Vysshaya shkola Publ., 1967, 600 p.

[17] Valtsiferov Yu.V., Dronov V.P. Chislennoe modelirovanie konvektivnogo tep-loobmena v tonkostennom cilindricheskom sosude s polusfericheskimi dnishchami, polnost'yu zapolnennom zhidkost'yu [Numerical simulation of con-vective heat transfer in a thin-walled cylindrical vessel with hemispherical bottoms completely filled with liquid]. Izvestiya Akademii nauk SSSR. Mekhanika zhidkosti i gaza [Proceedings of the USSR Academy of Sciences. Fluid and Gas Mechanics], 1984, no. 5, pp.204 -207.

[18] Verkin B. I., Kirichenko Yu. A., Rusanov K. V. Teploobmen pri kipenii krio-gennyh zhidkostej [Heat transfer during the boiling of cryogenic liquids]. Kiev, Nauk. Dumka Publ., 1987, 264 p.

[19] Patankar S. Numerical heat transfer andfluidflow. Washington, DC, Hemisphere Publishing Corporation, 1980, 210 p.

[20] Quazzani J., Garrabos Y. A new numerical algorithm for low Mach number supercritical fluid. Preprint Elsevier, 2007, 10 p.

[21] Samarskij A.A. Teoriya raznostnyh skhem [Theory of difference schemes]. Moscow, Nauka Publ., 1989, 616 p.

[22] Chung T.J. Computational fluid dynamics. Cambridge, Cambridge university press, 2010, 1012 p.

[23] Seo M., Jeong S. Analysis of self-pressurization phenomenon of cryogenic fluid storage tank with thermal diffusion model. Cryogenics, 2010, vol. 50, no. 9, pp. 549-555.

[24] Belyayev A.Yu., Ivanov A.V., Egorov S.D., Voyteshonok V.S., Mironov V.M. Pathways to solve the problem of cryogenic rocket propellant long storage in space. Proceedings of the International Aerospace Congress, Moscow, Russia, 1994, vol. 1, pp. 558-562.

[25] Grigoriev I.S., Meilikhov E.Z. Fizicheskie velichiny [Physical quantities]. Moscow, Energoatomizdat Publ., 1991, 1232 p.

[26] Solntsev Yu.P., Ermakov B.S., Sleptsov O.I. Materialy dlya nizkih i kriogennyh temperatur: enciklopedicheskij spravochnik [Materials for low and cryogenic temperatures: an encyclopedic reference]. Saint-Petersburg, Himizdat Publ., 2008, 768 p.

[27] Vargaftik N.B. Spravochnik po teplofizicheskim svojstvam gazov i zhidkostej [Handbook of thermophysical properties of gases and liquids]. Moscow, Nauka Publ., 1972, 720 p.

Gorodnov A.O., Research Associate of SSC Keldysh Research Center, research interests: fluid mechanics, thermophysics, numerical methods, mathematical modeling, cryogenic engineering. e-mail: an.ol.gorodnov@gmail.com