CC BY
46Математика к Математическое
моделирование
ХДК 536.2
Ссылка на статью:
// Математика и математическое моделирование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. №6. С. 44-60.
Б01:10.7463/шаШш.0615.0829350
Представлена в редакцию: 06.12.2015 © МГТУ им. Н.Э. Баумана
Математическое моделирование температурного состояния оболочки цилиндрической криогенной емкости при заполнении и опорожнении
Зарубин В. С.1'*, Зимин В. Н.1, Кувыркин Г. Н.1
1МГТУ им. Н.Э. Баумана, Москва, Россия
С учетом особенностей теплообмена в криогенной емкости построена математическая модель, описывающая температурное состояние тонкостенной оболочки цилиндрического бака для сжиженного газа при его заполнении и опорожнении. Проведен количественный анализ этой модели при постоянной скорости перемещения уровня криогенной жидкости, а также для случаев неподвижного уровня и уровня, совершающего гармонические колебания. Установлены предельные варианты квазистационарного распределения температуры вдоль образующей оболочки в подвижной системе координат при возрастании скорости заполнения или опорожнения емкости. Методом интегрального преобразования Лапласа решена нестационарная задача теплопроводности в подвижной системе координат для несмоченной части оболочки емкости. Анализ решения этой задачи позволил для случая перемещения уровня жидкости с постоянной скоростью оценить время установления квазистационарного распределения температуры вдоль образующей оболочки емкости. Представленные расчетные зависимости могут быть использованы при анализе напряженно-деформированного состояния и устойчивости оболочки криогенной емкости и для оценки потерь криогенной жидкости вследствие ее испарения.
Ключевые слова: криогенная емкость; модель температурного состояния; квазистационарное распределение температуры
Введение
Жидкие кислород и водород находят применение в качестве окислителя и горючего для жидкостных ракетнык двигателей [1, 2]. Сжиженный природный газ, основой которого является метан, рассматривается как перспективное моторное топливо для двигателей внутреннего сгорания [3, 4]. Одной из технических проблем, возникающих при использовании указанных криогенных жидкостей, является создание емкостей для их хранения, транспортировки и применения в двигательных установках. При проектировании и эксплуатации таких емкостей необходимо располагать достоверной информацией об их температурном состоянии, от которого зависят потери криогенных жидкостей за счет испарения и напряженно-деформированное состояние элементов конструкции емкостей.
Неравномерное распределение температуры вдоль образующей тонкостенной оболочки цилиндрического криогенного ракетного бака, локализованное в зоне уровня криогенной жидкости, приводит к искривлению оболочки и снижению допустимой осевой нагрузки, к возникновению опасности потери устойчивости оболочки при подготовке ракеты к старту и в полете с растущим ускорением [5]. Перемещение уровня криогенной жидкости в процессе заполнения или опорожнения бака при определенном сочетании параметров приводит к увеличению локальной неравномерности распределения температуры.
Наряду с экспериментальным изучением температурного состояния оболочки криогенной емкости необходимая при проектировании и отработке конструкции криогенных баков информация может быть получена методами математического моделирования [6, 7]. В данной работе с учетом особенностей теплообмена в криогенной емкости, в том числе при кипении криогенной жидкости на внутренней поверхности этой емкости, построена математическая модель, описывающая температурное состояние тонкостенной оболочки цилиндрического криогенного бака при его заполнении и опорожнении и представлен количественный анализ этой модели для случаев неподвижного уровня жидкости, его перемещения с постоянной скоростью и гармонических колебаний относительно некоторого среднего положения. Количественный анализ этой модели позволил установить предельные варианты квазистационарного распределения температуры вдоль образующей оболочки в подвижной системе координат при возрастании скорости заполнения или опорожнения емкости. Решение методом интегрального преобразования Лапласа нестационарной задачи теплопроводности в подвижной системе координат для несмоченной части оболочки емкости использовано для оценки времени установления квазистационарного распределения температуры в этой части оболочки.
1. Условия теплообмена в емкости
Интенсивность теплообмена смоченной и несмоченной участков внутренней поверхности оболочки криогенной емкости существенно различна, что связано, главным образом, с существенным различием плотности криогенной жидкости и газа, заполняющей верхнюю часть емкости. При неподвижном уровне жидкости процесс теплообмена на несмоченной части поверхности определяется естественной конвекцией газа. В случае сравнительно высокой скорости заполнения емкости при открытом дренажном клапане некторое влияние на интенсивность теплообмена может оказывать направленное движение газа, образующегося при испарении жидкости.
В качестве основной расчетной зависимости для определения среднего значения коэффициента теплообмена а на несмоченной вертикальной поверхности оболочки емкости можно использовать известную критериальную формулу [8]
Ш = С' Ба"', (1)
где Ки = а1/Хд — число Нуссельта; I — характерный размер поверхности, равный в данном случае длине образующей на несмоченном участке оболочки; Ад — коэффициент теплопроводности газа; Ба = Ог Рг, Ог = д0/ЗдI3 АТ/^2 и Рг = ид/ад — числа Рэлея, Грасгофа и Прандтля соответственно; д0 — ускорение поля тяготения (для неподвижной емкости на поверхности Земли д0 ~ 9,81 ); вд — коэффициент объемного расширения газа (в случае совершенного газа равен обратному значению его температуры Тд); АТ = |Тд — Тш|, Тад — температура поверхности теплообмена; ид и ад — соответственно кинематический коэффициент вязкости и коэффициент температуропроводности газа (значения всех тепло-физических характеристик газа следует принимать при температуре (Тд +Тад) /2). Параметры С' и и' зависят от числа Рэлея.
Условия эксплуатации криогенных емкостей таковы, что давление газа над уровнем жидкости обычно близко к атмосферному, а температура газа не превышает 273 К. В таких условиях при I « 1 м, АТ « 10 К и теплофизических характеристиках газов (азота, кислорода, метана, водорода, гелия) Ба > 2 ■ 107, что соответствует в формуле (1) значениям С' = 0,135 и и' = 1/3 [8], т.е. коэффициент а не зависит от характерного размера поверхности теплообмена. При этих условиях коэффициент теплопроводности азота и кислорода
Вт Вт
не превышает 0,025 ——, а метана — 0,032 —— [9], т.е. из формулы (1) следует оценка
„ м • К м К
Вт
а
>
(АТ)1/3 —2——, если значение АТ подставить в кельвинах. Для гелия и водорода
м2 ■ К
Вт Вт
Ад « 0,15 ... 0,18 [9] и а > 5(АТ)1/3 .
м ■ К ~ м2 ■ К
На смоченной части поверхности при неподвижном уровне жидкости с температурой Tf, близкой к ее температуре насыщения Т3 при существующем в емкости давлении, и близком к нулю значении АТШ = Тш — Tf процесс теплообмена также связан с механизмом естественной конвекции, а для оценки коэффициента теплообмена можно тоже использовать формулу (1) при значениях С' = 0,135 и и' = 1/3 и теплофизических характеристиках соответствующей криогенной жидкости. Для жидкого кислорода и метана коэффициент теплопровод-
Вт Вт
ности Аf при температуре насыщения равен соответственно 2,45 —— и 3,04 —— [10], что
м К м К
Вт Вт
приводит к оценке а > (90 ... 110)(АТад)1/3 . Для жидкого водорода Аf = 16,8 —— и
~ м2 - К м • К
а > 615(АТад)1/3 .
м2 ■ К
Из приведенных оценок следует, что в условиях естественной конвекции коэффициенты теплообмена на смоченной и несмоченной поверхностях отличаются примерно на два порядка, но пропорциональны одинаковой степени значений АТад и АТ. Близкие к формуле (1) оценки дает критериальная зависимость Ки = 0,0674(0г Рг1,29)1/3 [11], справедливая в интервалах 2,4 < Рг < 118 и 4 ■ 1010 < Ба < 9 ■ 1011 и полученная обработкой экспериментальных данных по теплообмену капельных жидкостей на вертикальных поверхностях при естественной конвекции.
По мере роста значения АТад возникает кипение криогенной жидкости (пузырьковый, переходной и пленочный режимы [12, 13]). В случае так называемого недогрева жидкости,
когда Tf < Т8, т.е. в емкости находится переохлажденная криогенная жидкость, кипение возникает при более высоких значениях АТт [11]. Такое же влияние может оказывать и движение уровня жидкости при заполнении или опорожнении емкости.
Для количественного описания интенсивности теплообмена при кипении предложено достаточно много математических моделей [11, 12] и выполнено большое число экспериментальных исследований, обработка результатов которых представлена эмпирическими зависимостями, включающими значительное количество параметров. Важными характеристиками процесса кипения являются ординаты экстремумов на кривой зависимости плотности q теплового потока от АТШ. На рис. 1 в качестве примера представлен в логари-фимических координатах график такой зависимости при кипении воды [10]. Пунктирная линия соответствует режиму естественной конвекции без кипения, сплошная кривая между точками А и В — режиму пузырькового кипения, штриховая между В и С — переходному, а штрихпунктирная — пленочному режиму кипения. Переходный режим является неустойчивым. В точке В последующее увеличении АТт вызывает переход пузырькового режима к пленочному Наоборот, уменьшение АТт в точке С приводит к переходу от пленочного режима кипения к пузырьковому.
Рис. 1. Зависимость плотности теплового потока от разности температур
Ординату максимума зависимости q от АТт, соответствующего первому кризису кипения, можно вычислить по формуле [10, 11]
q* = (0,13т(адо- Рд)р2д)1/4 + 0,0012™^рд)1/д) (1 + 0,065^/рд)0,8^(Тв - Tf )/т),
где т, а, pf, Cf — соответственно теплота испарения, поверхностное натяжение, плотность и удельная массовая теплоемкость жидкости, рд и V — плотность газа и скорость перемещения уровня жидкости в емкости. Для оценки ординаты минимума зависимости q от АТШ, соответствующего второму кризису кипения, рекомендовано соотношение [11]
q* = (0,11... 0,14)трд(аgo(pf - рд)/р))
1/4
для оценки абсциссы этого минимума в случае кипения на вертикальной поверхности [12] —
АТ* = * = (0,083... 0,089) ^ (^ — рд )рд ^)1/3 ( " )1/2. а* Ад V (рf + рд)2 ) Vgo(Pf — рд))
Коэффициент теплообмена при пузырьковом кипении в условиях естественной конвекции можно оценить по формуле
ак = А, , /
г2,5
« 1п-4т-, -0,35 / Л ,гт \ 0,Л 10/3
7 ■ 10 4 Рг- , /pАf АТ^ ,
(^5д0^ — Рд))0,Ч грдаf
полученной преобразованием известной формулы Кутателадзе [11, 14]. Здесь Аf и Pгf = = ^/af — коэффициент теплопроводности и число Прандтля жидкости; ^ и af — соответственно кинематический коэффициент вязкости и коэффициент температуропроводности
жидкости; р — давление. При кипении водорода и кислорода экспериментально измерены
Вт Вт
значения ак соответственно до 5 ■ 105 -и 2 ■ 105 -. В случае пленочного режима ки-
м2 ■ К м2 ■ К
пения на вертикальной поверхности для оценки коэффициента теплообмена рекомендовано
( ) 1/3
соотношение [12] апл = 0,25Ад(#0^/рд — 1)/(ад^ .
2. Построение математической модели
Условия теплообмена цилиндрической оболочки вертикально расположенной цилиндрической криогенной емкости в зоне перехода через уровень жидкости претерпевают резкое изменение, что приводит к возникновению существенно неравномерного распределения температуры Т вдоль образующей оболочки в этой зоне. Изменение этого распределения во времени £ и вдоль образующей оболочки по координате г описывает дифференциальное уравнение с частными производными [15]
зт(м ьдЛдОМ
^ + а(г,£)(Т *(г,£) — Т (г,£)) + а*(Т* — Т (*,£)), (2)
где с, р и А — удельная массовая теплоемкость, плотность и коэффициент теплопроводности материала оболочки с постоянной толщиной Л; а — зависящий (в общем случае перемещения уровня жидкости) от координаты и времени коэффициент теплообмена на внутренней поверхности оболочки со средой, имеющей температуру Т*, также в общем случае являющейся функцией координаты и времени; а* — коэффициент теплообмена на внешней поверхности оболочки со средой, имеющей температуру Т*. Положительное направление координаты г выберем вертикально вверх.
Для однозначного решения уравнения (2) помимо начального распределения температуры Т(г, 0) в момент времени £ = 0, принимаемый за начальный, необходимо задать граничные условия. Если условия теплообмена на смоченной и несмоченной поверхностях оболочки определяют постоянные значения соответственно а1, Т* и а2, Т2*, то по мере удаления от уровня жидкости роль передачи теплоты вдоль образующей оболочки путем
теплопроводности становится несущественной [ 16] (формально дТ/дг ^ 0 при г ^ -то и дТ/дг ^ 0 при г ^ то), а температура оболочки стремится к значениям
гр( П а1Т1* + а*Т* , П а2Т2 + а*Т*
Т (г,с) ^ Т1 =-, г ^ —то; Т (г,с) ^ Т2 =-, г ^ —то.
а1 + а* ад + а*
Тогда граничные условия можно задать в виде
Т(г,С) = Т^ г ^—то; Т(г,С) = Тд, г ^ то. (3)
В случае неподвижного уровня жидкости, положение которого соответствует значению г = 0, и независящего от температуры коэффициента теплопроводности Л при постоянных во времени условиях теплообмена установившиеся распределения температуры Т1 (г) и Т2 (г) в смоченной и несмоченной частях оболочки будет описывать вытекающая из равенства (2) система двух обыкновенных дифференциальных уравнений второго порядка
^ + а2(Т1 - ВД) =0, г е (-то, 0), (4)
+ ад2(:Гд - Тд(г)) =0, г е (0, то), (5)
где ад = (а1 + а*)Л,/Л и ад = (ад + а*)Л,/Л. Решение системы уравнений (4) и (5) при граничных условиях (3) и условиях
-T
Ti(0) = T2(0)= To, —1
dz
z=0
-T2
dz
(6)
z=0
непрерывности при г = 0 температуры и теплового потока имеет вид [5, 15]
Ti(z) - T = exp(o;iz/h) z < о. T2 - T2(z) = exp(-a2z/h) z > 0 (7)
T2 - Tx = 1 + ai/02 , z . T2 - Tx = 1 + Ö2/Ö1 , z ' ()
Равенства (7) могут быть использованы при формулировке начального условия для уравнения (2), если началу перемещения уровня жидкости предшествует достаточно длительный период неподвижного положения этого уровня.
Из равенств (7) следует формула для температуры T0 при z = 0 на границе между смоченной и несмоченной частями оболочки. В безразмерном виде эту температуру можно представить в виде
To - Tl = а ' (8) T2 - Тх ai + a2
Поскольку exp(-3) ~ 0,05, практически можно считать, что неравномерность установившегося распределения температуры в оболочке достаточно учитывать лишь при -3/ах < < z/h < 3/a2. Коэффициенты теплообмена a2 и а* обычно сопоставимы по значениям, а отношение ах/а2 достаточно велико в силу высокой интенсивности теплообмена на смоченной поверхности оболочки. Это приводит к тому, что значение температуры T0 мало отличается от Тх. В таком случае температуру всей смоченной части оболочки в первом
приближении можно считать равной температуре жидкости, т.е. положить Т1(г) ~ Т1 ~ Т*, а распределение температуры в несмоченной части оболочки представить в виде
Т2(г) « Т2 — (ТТ2 — Т*) ехр^ — (9)
Такое допущение существенно упрощает количественный анализ установившегося распределения температуры в оболочке.
3. Переход к подвижной точке отсчета координаты
При изменении условий теплообмена в оболочке возникает нестационарный процесс теплопроводности, описываемый уравнением (2). Одной из причин локального изменения условий теплообмена является движение уровня жидкости при заполнении или опорожнении емкости. В случае постоянной скорости V перемещения этого уровня целесообразно в уравнении (2) перейти к координате ( = г — V*, начало отсчета которой связано с текущим положением уровня жидкости (V > 0, если происходит заполнение емкости). В результате вместо уравнения (2) получим систему двух уравнений с частными производными
1 та + + (%у2(Т,—тЫ. с<0, (10)
а д* 2 а д( \ Л
1 ОД«,*) д2Т2(С,£) , V дТ^,*) , (а2'
2
а в* *> + а^р+ ^ ^—Т2(И, с>(11)
Граничные и начальные условия для уравнений (10) и (11) будут такими же, как и для уравнения (2), но с учетом замены координаты г координатой (, а именно, вместо условий (3) получим
т «,*) = 7, с т «,*) = Т2, с (12)
вместо (6) —
Т1(0,*) = Т2(0,*),
= ¿7
с=0 ¿С
, (13)
С=0
а для задания начального распределения температуры можно использовать формулы (7) с указанной заменой координаты. В этом состоит одно из преимуществ перехода к координате (. Второе преимущество заключается в том, что в случае граничных условий, не зависящих от * в подвижной системе координат, 5Т1((,*)/5* ^ 0 и дТ2(£,*)/$* ^ 0 при * ^ то. Это означает, что с течением времени после начала движения уровня жидкости нестационарные распределения Т1 ((, *) и Т2 ((, *) температуры в оболочке будут стремиться к некоторым установившимся, определяемым из решения системы двух обыкновенных дифференциальных уравнений
2
2 ' а V Л /
^ + аТ + (IГ(72 — Т2(С)) О °. (.5)
которые следуют из уравнений (10) и (11), если их правые части приравнять нулю.
^++(!)>■—т^ с<°. (и)
Система уравнений (14), (15) описывает так называемое квазистационарное распределение температуры в оболочке емкости, не зависящее от времени в подвижной системе координат Решение этой системы с использованием граничных условий (12)и(13) совпадает по структуре с соотношениями (7) и (8), но при условии замены параметра а1 на
При сравнительно большой скорости заполнения емкости, когда 2аа2/(уЬ) ^ 1, роль передачи теплоты вдоль образующей оболочки путем теплопроводности мала. Тогда для смоченной части оболочки [5]
а для несмоченной — Т2(() ~ Т2. Наоборот, при быстром опорожнении емкости, когда 2аа1 /(уЬ) ^ 1, для смоченной части оболочки Т1 (С) ~ Т, а для несмоченной —
Эти соотношения устанавливают крайние возможные квазистационарные распределения температуры вдоль образующей оболочки. На рис. 2 приведено взаимное расположение точек отсчета координат £ и £ в случае опорожнения емкости (у < 0) и показан качественный характер распределений температуры вдоль образующей оболочки при у < 0 и у > 0, а также в предельных случаях у ^ —го и у ^ го.
и параметра а2 на
Рис. 2. Распределения температуры в оболочке при движении уровня жидкости
4. Решение нестационарных задач
Помимо поступательного перемещения уровня жидкости в емкости возможно колебание этого уровня относительно некоторого среднего неподвижного положения. При достаточно высокой частоте колебаний можно не учитывать передачу теплоты теплопроводностью вдоль образующей оболочки и, возвращаясь к неподвижной системе координат, в случае фиксированных амплитуды z0 и круговой частоты и колебаний периодическое изменение температуры Т(z,t) во времени в фиксированном сечении с координатой z е (—z0,z0) будет описывать обыкновенное дифференциальное уравнение первого порядка
cphdTdr) = a(t)(f(t) - T(z,t)), (16)
где a(t) и T(t) — изменяющиеся скачкообразно по времени приведенные коэффициент теплообмена и температура среды. Если — z0 ^ z ^ z0 sin ut, то в данный момент времени a(t) = a1 + a* и T(t) = Ть а в противном случае a(t) = a2 + a* и T(t) = Т2. Следовательно, в течение периода t* = 2n/u колебаний условий теплообмена для фиксированного значения z в промежутке времени t1 = t* arcsin z1, где z1 = (z + z0)/(2z0) уравнение (16) имеет вид
CphdTdM = (ai + a*) (Tfi — t (z, t)), (17)
а в промежутке времени t2 = t* — t1 —
= (a2 + a*) (T2 — T (z, t)). (18)
Начальной для уравнения (17) будет температура T2 в конце промежутка времени t2, определяемая из решения уравнения (18) с начальной температурой T1 < T2 в конце промежутка времени t1, определяемой из решения уравнения (17). Эти температуры можно найти из системы уравнений
Т/ р гр rpf
1 — Т ( \ T 2 — Т 2 / \
гр, р = exp( —^ Тр-= eXP( — Т2 ),
Т2 — Т1 Т 2 — Т1
(«1 + a*)Í1 («2 + a*)Í2 ~ ,
где T1 =---и n =---. Отсюда следуют расчетные формулы
О T2 — Т1 1 — exp( — т2) T1 — Т1 q , \
= —-— =--------, —-— = exp( — ti).
T2 — Т1 1 — exp(—n) exp(—T2) T2 — T1
После вычисления предельных значений температур Т, и Т2 можно представить решение уравнений (17) и (18) в виде
Т1 — Т / T1t \ п Т2 — Т2 / T2t \
—-^ = exp —- , t е [0, t1]; -J-2 = exp —- , t е [0, t2].
T2 — T1 V t^' L' 1J' T2 — T1 V t^ L' 2J
Вернемся к случаю поступательного перемещения уровня жидкости с постоянной скоростью v. Для оценки периода времени, по истечении которого нестационарное распределение
температуры в оболочке с определенной точностью можно считать совпадающим с квазистационарным распределением, рассмотрим упрощенную задачу, приняв в силу высокой интенсивности теплообмена оболочки с жидкостью Т2(0, ¿) = 7\(0, ¿) = 7\. В этом случае следует решить уравнение (11) с граничными условиями Т2(0, ¿) = Т\ и Т2(£2,¿) ^ Т2 при ( ^ то. В качестве начального условия примем распределение температуры, определяемое равенством
Т2 - Т2(С) А а2С А
ехр(^ —^ , г ^ 0. (19)
7 - П
Для решения этой задачи применим интегральное преобразование Лапласа [17,18,19]. Тогда уравнение (11) в сочетании с начальным условием примет вид
+ V - (в + ( 02 УН = 7 - 7 /ос А - (£ + ( а 72 (20)
2 а у а V Н / у а V Н / \а Ц / / « '
где так называемое изображение искомой функции (оригинала) Т2((, ¿) равно
72«) = / Т2«,*)ехр(-в*) ^
0
Граничными условиями для полученного относительно изображения обыкновенного дифференциального уравнения (20) будут Т2(0) = Т^/в и Т2(() ^ Т2/в при ( ^ то. Из второго граничного условия следует, что решение уравнения (20) должно быть ограниченным при ( ^ то, т.е. из двух корней характеристического уравнения, соответствующего однородному уравнению, сохраняет смысл лишь корень
V
г = — ---
2а \
V2 /«2\2 в
(2а)2 + V Н) + а.
Тогда решение неоднородного уравнения (20) можно представить в виде
Ф п (г А . Т Т2 - А а2С
Т 2 = д ехр — +------ ехр
в в + ^2/Н V Н Из граничного условия при £ = 0 получим
Т2 - Т 1 Тт2 - 7
в + ^2/Н в В итоге решение задачи в изображениях примет вид
т2(С) = ( 72 - 7 - Тт2 - 7 \ ехр(^л + 72 - ттъ - 7 / оС\ (21)
у в + ^2/Н в ) в в + va2/Н V Н /
Из равенства (21), используя таблицы соответствия изображений и оригиналов [18], получаем
7Ж, *) = 7 - (7 - Т1) ехр [- ^ С^ + (7 - ((, *) ехр [-2^) , (22)
где
^(С,*) = ^(С,*)ехр(—ис^/Ь) - ¿2«,*)
*Ж,*) = 1 ехр (— ^ + М) егГе (- (V — 02е-^ + ПЧ' ; 2 ^ 2а Ь ] \2Vat \2а Ь) )
+ 1ехр( < — М1 егГе(^ + (^ — %) 2 Р\2а Ь ) \2а Ь) )
ЯК, О = 2ехр(—Уё^/Й^с/2) ^ ^ — ЛР!
+ 2 ехр У (и/а)2 + (2а2/Ь)2(/^ eгfc (
1 Г,—ГТ^-^-/<Л ^ С /и2* а|а*
1 —^^л/а^ " V 4а ' Ь2
Если функция ^(С,*) ограничена при £ ^ то, то полученное решение удовлетворяет граничному условию Т2(£, *) ^ Т2 при £ ^ то.
По мере возрастания * второе слагаемое в правой части равенства (22), характеризующее
начальное распределение температуры, убывает пропорционально множителю ехр I--— I.
При этом должна возрастать роль третьего слагаемого в правой части этого равенства, определяемая свойствами функции ^(С,*). На эти свойства существенное влияние оказывает функция егТс(ж), дополняющая до единицы функцию ошибок егТ(ж) Гаусса и ограниченная в интервале (—то, то), причем
2 ^
егБс(ж) = 1 — егДж) = —^ J ехр(—£2) егБс(то) = 0, ег4С(0) = 1, егБс(-то) = 2.
Отсюда, в частности, следует, что ^(£, 0) = 0, т.е. полученное решение удовлетворяет начальному условию.
В силу ограниченности функции егТс(ж) функция *) ехр(—иа2*/Ь) убывает по мере возрастания * в силу убывания экспоненциального множителя. При * ^ то второе слагаемое в правой части равенства (23) обращается в нуль, а первое слагаемое принимает вид
с (у\2 2
ехР —пУ " +
2 а Ь
т.е. полученное решение (22) в этом случае описывает в подвижной системе координат установившееся (квазистационарное) распределение температуры
7-2(0 = Т2 — й — Т1) ехр^—2 (а) — Ш2 +Ш2). (24)
Если использовать принятые при решении уравнения (11) граничные условия, то такое распределение температуры непосредственно следует из решения обыкновенного дифференциального уравнения в виде равенства нулю правой части уравнения (11).
Для количественного анализа полученное решение (22) представим в виде
©(V, С, С) = Т ,С) = ехр(-( - С) - Р(V , С, С) ехр(-*С):
Т2- Т
где
V =
2аа2'
с =
а2 а£
С =
«2 С
к2 к2 к
Р(V, С С) = Р(V,С С) ехр(-2^") - Р2(V,С С);
Р (V, С, С) = ^ ехр( (1 - V)С) егГсГ= + (1 - V)-С +
+ ^ ехр( -(1 - ег^- (1 - V)-
Р2(V, С, С) = ^
ехр^-—1 + V2С^ ег!с^С
-1 + V2 VС | +
1
+ ^ехр
(—1 + V2 С) ег!"с^ с + —1 + V2 VС).
На рис. 3 в полулогарифмической системе координат приведены результаты проведенного количественного анализа в виде графиков распределений относительной безразмерной температуры В вдоль безразмерной координаты £ с подвижной точкой отсчета, соответствующей перемещающемуся уровню жидкости. Начальное распределение температуры представлено сплошной кривой со светлыми кружками. Скорость перемещения уровня входит в безразмерный параметр V, для нескольких значений которого на рисунке приведены квазистационарные распределения температуры, рассчитанные по формуле (24) (кривая со светлыми треугольниками построена при V = 100, со светлыми ромбами и квадратами — соответственно при V = 20 и V = 5, с темными кружками — при V = 1, с темными трегольниками и ромбами — соответственно при V = -1 и V = -2, т.е. для случая опорожнения емкости в отличие от положительных значений V, относящихся к случаю заполнения емкости).
Тонкая сплошная кривая около каждой из кривых с символами описывает на рис. 3 нестационарное распределение температуры для фиксированного значения 7, подобранного с
Рис. 3. Результаты количественного анализа
использованием решения (22) из условия практического совпадения нестационарного распределения температуры с соответствующим квазистационарным. Каждое из значений ¿* можно рассматривать как оценку периода времени (в безразмерном виде), необходимого для завершения процесса перестройки распределения температуры вдоль образующей оболочки от начального до квазистационарного. Установлено, что по мере убывания параметра V соответствующие им значения ¿* монотонно возрастают: ¿* = 0, 02 при V = 100, i* = 0,1 при V = 20, i* = 0, 35 при V = 5, i* = 1 при V =1. При двух отрицательных значениях V обе тонкие сплошные кривые построены для одинакового значения ¿* = 1, 8. При этом заметное отклонение от соответствующих квазистационарных распределений выявлено лишь при в > 0,8 в зоне, достаточно удаленной от движущегося уровня жидкости.
Заключение
Анализ математической модели, построенной с учетом особенностей теплообмена на внутренней поверхности криогенной емкости (включая особенности различных режимов кипения криогенной жидкости на этой поверхности) позволил получить расчетные зависимости для предельных вариантов распределения температуры вдоль образующей цилиндрической оболочки емкости при неподвижном уровне жидкости, его поступательном движении с постоянной скоростью и в случае гармонических колебаний относительно некоторого среднего положения. Решение нестационарной задачи теплопроводности в несмоченной части оболочки использовано для оценки времени установления квазистационарного распределения температуры в этой части оболочки.
Работа выполнена по гранту НШ-1432.2014.8 программы Президента РФ государственной поддержки ведущих научных школ и в рамках проекта 1712 в сфере научной деятельности в части государственного задания №2014/104 Минобрнауки РФ, а также в рамках государственного задания по проекту № 1.2640.2014.
Список литературы
1. Феодосьев В.И. Основы техники ракетного полета. М.: Наука, 1979. 496 с.
2. Ковалев Б.К. Развитие ракетно-космических систем выведения. М.: Изд-во МГТУ им.
Н.Э.Баумана, 2014. 400 с.
3. Альтернативные топлива для двигателей внутреннего сгорания / Под общ. ред.
А.А. Александрова, В.А. Маркова. М.: ООО НИЦ "Инженер", ООО "Онико-М", 2012.
791 с.
4. Горбачев С.П., Коледова К.И., Красноносова С.Д. Термодинамические модели заправки
резервуара криогенной жидкостью // Технические газы. 2011. № 5. C. 32-40.
5. Балабух Л.И., Колесников К.С., Зарубин В.С., Алфутов Н.А., Усюкин В.И., Чижов В.Ф.
Основы строительной механики ракет. М.: Высшая школа, 1969. 496 с.
6. Зарубин B.C., Кувыркин Г.Н. Математическое моделирование термомеханических процессов при интенсивном тепловом воздействиии // Теплофизика высоких температур. 2003. Т. 41, № 2. C. 300-309.
7. Зарубин B.C., Кувыркин Г.Н. Особенности математического моделирования технических устройств//Математическое моделирование и численные методы. 2014. Т. 1, № 1-1. C. 517.
8. Кутателадзе C.C., Боришанский В.М. Справочник по теплопередаче. М.: Госэнергоиздат, 1958. 414 с.
9. Физические величины: ^равочник / Под ред. H.C. Григорьева, Е.З. Мейлихова. М.: Энергоатомиздат, 1991. 1232 с.
10. Малков М.П., Данилов И.Б., Зельдович А.Г., Фрадков А.Б. ^равочник по физико-техническим основам криогеники / Под ред. М.П. Малкова. М.: Энергоатомиздат, 1985. 432 с.
11. Теория тепломассообмена/Под ред. А.И. Леонтьева. М.: Изд-воМГТУ им. Н.Э. Баумана, 1997. 683 с.
12. Григорьев В.А., Павлов Ю.М., Аметистов Е.В. Кипение криогенных жидкостей / Под ред. Д А. Лабунцова. М.: Энергия, 1977. 289 с.
13. Веркин Б.И., Кириченко Ю.А., Русанов К.В. Теплообмен при кипении криогенных жидкостей. Киев: Наукова думка, 1987. 240 с.
14. Кутателадзе C.C. Основы теории теплообмена. Новосибирск: Наука, 1970. 659 с.
15. Зарубин В^. Температурные поля в конструкции летательных аппаратов. М.: Машиностроение, 1966. 216 с.
16. Черкасов C.E, Моисеева Л.А. Влияние продольного перетока тепла на распределение температуры в движущемся ребре при скачкообразном распределении коэффициента теплообмена // Теплофизика высоких температур. 2015. Т. 53, № 5. C. 807-809.
17. Карслоу Г., Егер Д. Теплопроводность твердых тел. М.: Наука, 1964. 488 с.
18. Лыков А.В. Теория теплопроводности. М.: Высшая школа, 1967. 600 с.
19. Карташов Э.М. Аналитические методы в теории теплопроводности твердых тел. М.: Высшая школа, 2001. 550 с.
Mathematics i Mathematical Modelling
Electronic journal of the Bauman MSTU
Mathematics and Mathematical Modelling of the Bauman MSTU, 2015, no. 6, pp. 44-60.
DOI: 10.7463/mathm.0615.0829350
Received: 06.12.2015
© Bauman Moscow State Technical University
http://mathmjournal.ru
Mathematical Modeling of the Thermal Shell State of the Cylindrical Cryogenic Tank During Filling and Emptying
Zarubin V.S.1*, Zimin V.N.1, Kuvyrkin G.N.]
1 Bauman Moscow State Technical University, Russia
Keywords: cryogenic tank, model of the temperature state, quasi-stationary temperature distribution
Liquid hydrogen and oxygen are used as the oxidizer and fuel for liquid rocket engines. Liquefied natural gas, which is based on methane, is seen as a promising motor fuel for internal combustion engines. One of the technical problems arising from the use of said cryogenic liquid is to provide containers for storage, transport and use in the propulsion system. In the design and operation of such vessels it is necessary to have reliable information about their temperature condition, on which depend the loss of cryogenic fluids due to evaporation and the stress-strain state of the structural elements of the containers.
Uneven temperature distribution along the generatrix of the cylindrical thin-walled shell of rocket cryogenic tanks, in a localized zone of cryogenic liquid level leads to a curvature of the shell and reduce the permissible axle load in a hazard shell buckling in the preparation for the start of the missile in flight with an increasing acceleration. Moving the level of the cryogenic liquid during filling or emptying the tank at a certain combination of parameters results in an increase of the local temperature distribution nonuniformity.
Along with experimental study of the shell temperature state of the cryogenic container, methods of mathematical modeling allow to have information needed for designing and testing the construction of cryogenic tanks. In this study a mathematical model is built taking into account features of heat transfer in a cryogenic container, including the boiling cryogenic liquid in the inner surface of the container. This mathematical model describes the temperature state of the thin-walled shell of cylindrical cryogenic tank during filling and emptying. The work also presents a quantitative analysis of this model in case of fixed liquid level, its movement at a constant speed, and harmonic oscillations relative to a middle position. The quantitative analysis of this model has allowed to find the limit options for quasi-stationary temperature distribution along the shell generatrix in the moving coordinate system with an increase in the rate of filling or emptying
the tank. Solution of a non-stationary heat conduction problem in moving coordinate system for unwetted part of the shell containers by the integral Laplace transform method is used to estimate the time required to define a quasi-stationary temperature distribution in this part of the shell.
References
1. Feodos'ev V.I. Osnovy tekhniki raketnogo poleta [The bases of techniques of rocket flight]. Moscow, Nauka Publ., 1979, 496 p. (in Russian).
2. Kovalev B.K. Razvitie raketno-kosmicheskikh sistem vyvedeniya [The development of rocket and space systems launch]. Moscow, Bauman MSTU Publ., 2014, 400 p. (in Russian).
3. Aleksandrov A.A., Markov V.A., ed. Al'ternativnye topliva dlia dvigatelei vnutrennego sgo-raniia [Alternative fuels for internal combustion engines]. Moscow, Ltd SIC "Engineer", LLC "Oniko-M", 2012, 791 p. (in Russian).
4. Gorbachev S.P., Koledova K.I., Krasnonosova S.D. Thermodynamic Models of a Filling Tank by the Cryogenic Liquid. Tekhnicheskie gazy = Industrial Gases, 2011, no. 5, pp. 32-40. (in Russian).
5. Balabukh L.I., Kolesnikov K.S., Zarubin V.S., Alfutov N.A., Usyukin V.I., Chizhov V.F. Osnovy stroitel'noi mekhaniki raket [Fundamentals of building mechanics of rockets]. Moscow, Vysshaya shkola Publ., 1969, 496 p. (in Russian).
6. Zarubin V.S., Kuvyrkin G.N. Mathematical modeling of thermomechanical processes under intense heat exposure. Teplofizika vysokikh temperatur = High Temperature, 2003, vol.41, iss. 2, pp. 300-309. (in Russian).
7. Zarubin V.S., Kuvyrkin G.N. Features of mathematical modeling of technical devices. Matem-aticheskoe modelirovanie i chislennye metody = Mathematical modeling and numerical methods, 2014, vol. 1, no. 1. pp. 5-17. (in Russian).
8. Kutateladze S.S., Borishanskii V.M. Spravochnikpo teploperedache [Handbook of heat transfer], Moscow, GosenergoizdatPubl., 1958, 414 p. (in Russian).
9. Grigor'ev I.S., Meilikhov E.Z., ed. Fisicheskie velichiny [Phisical values]. Moscow, Energoat-omizdatPubl., 1991, 1232 p. (in Russian).
10. Malkov M.P., Danilov I.B., Zel'dovich A.G., Fradkov A.B. Spravochnik po fiziko-tekhnicheskim osnovam kriogeniki [Handbook of physical and technical fundamentals of cryogenics]. Moscow, EnergoatomizdatPubl., 1985. 432 p. (in Russian).
11. Isaev S.I., Kozhinov I.A., Kofanov V.I., Leont'ev A.I., eds. Teoriia teplomassoobmena [The theory of heat and mass transfer]. Moscow, Bauman MSTU Publ., 1997. 683 p. (in Russian).
12. Grigoriev V.A., Pavlov Yu.M., Ametistov E.V., Labuntsov D.A., eds. Kipenie kriogennykh zhidlostei [Boil cryogenic liquids]. Moscow, Energiia Publ., 1977, 289 p. (in Russian).
13. Verkin B.I., Kirichenko Yu.A., Rusanov K.V. Teploobmenpri kipenii kriogennykh zhidlostei [Boiling heat transfer of cryogenic liquids]. Kiev, Naukova Dumka Publ., 1987. 240 p. (in Russian).
14. Kutateladze S.S. Osnovy teorii teploobmena [The bases of the theory of heat transfer]. Novosibirsk, NaukaPubl., 1970, 659 p. (in Russian).
15. Zarubin V.S. Temperaturnyepolia vkonstruktsii letatel'nykh apparatov [The temperature field in the construction of aircraft]. Moscow, Engineering publ., 1966, 216 p. (in Russian).
16. Cherkasov S.G., Moiseyeva L.A. Influence of longitudinal flow heat on the temperature distribution in the moving edge with spasmodic distribution of the heat transfer coefficient. Teplofizika vysokikh temperatur = High Temperature, 2015, vol. 53, iss. 5, pp. 807-809. (in Russian).
17. CarslawH.S., Jaeger J.C. Conduction of Heat in Solids. 2nd ed. Oxford University Press, 1959. (Russ. Ed. Karslou G., Eger D. Teploprovodnost' tverdykh tel. Moscow, Nauka Publ., 1964. 488 p.).
18. Lykov A.V. Teoriia teploprovodnosti [The theory of heat conduction]. Moscow, Vysshaya shkola Publ., 1967. 600 p. (in Russian).
19. Kartashov E.M. Analiticheskie metody v teorii teploprovodnosti tverdykh tel [Analytical methods in the theory of heat conduction of solids]. Moscow, Vysshaya shkola Publ., 2001, 550 p. (in Russian).
