Квазиодномерная нестационарная модель процессов в морских газопроводах Текст научной статьи по специальности «Математика»

УДК 532.517+532.542

Н. Н. Ермолаева, Г. И. Курбатова

Вестник СПбГУ. Сер. 10. 2015. Вып. 3

КВАЗИОДНОМЕРНАЯ НЕСТАЦИОНАРНАЯ МОДЕЛЬ ПРОЦЕССОВ В МОРСКИХ ГАЗОПРОВОДАХ

Санкт-Петербургский государственный университет, Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7/9

Предложены математические модели установившихся и неустановившихся процессов в потоке газа и процессов теплообмена газа с окружающей средой. Модели включают учет возможного оледенения внешней поверхности газопровода в северных морях. Представлены расчеты по ряду упрощенных моделей. Приведены качественные и количественные оценки допустимости упрощенных моделей в случае медленных изменений во времени температуры газового потока. Получена оценка максимальной скорости нарастания льда для упрощенной модели процессов. Библиогр. 9 назв. Ил. 1.

Ключевые слова: газопроводы, образование льда, анализ термодинамической модели.

N. N. Ermolaeva, G. I. Kurbatova

QUASI-ONE-DIMENSIONAL NON-STATIONARY MODEL OF PROCESSES IN A SEA GAS-PIPELINES

St. Petersburg State University, 7/9, Universitetskaya embankment, St. Petersburg, 199034, Russian Federation

The article presents the mathematical model of non-stationary heat exchange processes between the gas flow in the sea gas-pipeline and the envirinment. The model takes into consideration the possibility of pipeline glaciation. The simplified variants of the model are discussed. They provide insight into the qualitative and the quantitative assessments of the admissibility of different simplifications. Refs 9. Fig. 1.

Keywords: gas-pipelines, ice formation, analysis of thermodynamic model.

Для описания моделей процессов в потоке газа и процессов теплоообмена с окружающей средой были приняты следующие обозначения: r, z — радиальная и осевая координаты в цилиндрической системе координат (r,p,z); t — время; rx, tx — характерные значения r, t; R — внутренний радиус газопровода; Si — толщина первого слоя обшивки газопровода; S2 — толщина второго слоя обшивки; Ri = R + Si; R2 = R + S1 + S2 — внешний радиус газопровода. Индексы: 1 — область первого слоя обшивки (стали); 2 — область второго слоя обшивки (бетона); l — область слоя льда на внешней поверхности газопровода; v — область эффективного теплового по-гранслоя воды; Xi — коэффициент теплопроводности в i-й области (i = 1, 2, l, v); а — коэффициент теплоотдачи на внутренней поверхности газопровода; pi — плотность материала i-й области (i = 1, 2, l, v); ci — удельный коэффициент теплоемкости материала i-й области (i = 1, 2, l, v); ai — безразмерный коэффициент температуропроводности в i-й области (i = 1, 2, l, v); Ti(r,t) — распределение температуры в i-й

Ермолаева Надежда Николаевна — кандидат физико-математических наук, доцент; e-mail: n.ermolaeva@spbu.ru

Курбатова Галина Ибрагимовна — доктор физико-математических наук, профессор; e-mail: gi_kurb@mail.ru

Ermolaeva Nadeczda Nikolaevna — candidate of physical and mathematical sciences, docent; e-mail: n.ermolaeva@spbu.ru

Kurbatova Galina Ibragimovna — doctor of physical and mathematical sciences, professor; e-mail: g.kurbatova@spbu.ru

области (i = 1, 2, l, v); T0(r) — начальное распределение температуры в i-й области (i = 1, 2, l, v); Si(t) — толщина слоя льда на внешней поверхности газопровода; ¿ю — начальная толщина слоя льда; p(z,t) — плотность газового потока; u(z,t) — скорость газового потока; p(z,t) — давление в газовом потоке; T(z,t) — температура газа; e(z,t) — массовая плотность полной энергии газа; e(z,t) — массовая плотность внутренней энергии газа; Л — коэффициент гидравлического сопротивления потока газа; g cos a(z) — проекция вектора ускорения силы тяжести на ось z; T* — температура фазового перехода вода-лед для заданной солености морской воды; T* — температура морской воды на удалении от газопровода; Q — удельная теплота плавления льда; S* — толщина эффективного теплового погранслоя воды; S*i — толщина эффективного теплового погранслоя воды при наличии слоя льда; L = ^^ (г£р) — оператор Лапласа в цилиндрической системе координат (г, z, ф) при = щ = 0; q(r,t) = Aj^- — поток тепла (радиальная составляющая вектора плотности потока внутренней энергии); qo(r) — поток тепла в установившихся процессах; Ai, Bi — константы в стационарной зависимости Ti(r) для i-й области (i = 1, 2, l, v); Ar, At — безразмерные величины шагов по радиусу и по времени в алгоритме численного интегрирования; Ri = R2 + Si (t) — внешний радиус газопровода с учетом оледенения; Dt — интервал времени; ASi — приращение слоя льда за Dt; t — момент времени понижения внешней поверхности газопровода до температуры T*; tk, Tk, T^ — параметры в законе поведения T(t); w(z,t) — объемная плотность в единицу времени распределенных источников (стоков) внутренней энергии в потоке газа.

Постановка задачи. Для протяженных морских газопроводов, работающих без промежуточных подстанций, характерны большие расходы и сверхвысокие давления. Поэтому газ в части газопровода может иметь температуру ниже температуры фазового перехода T* вода-лед. Для южных морей это не приводит к оледенению газопровода. Иначе обстоит дело в северных морях, где нередко температура воды близка к температуре T*, потому при остывании потока ниже определенной температуры T (T < T*) оледенение газопровода становится возможным. Целью настоящей работы является исследование теплообмена между потоком газа и окружающей средой для морских газопроводов, эксплуатируемых в условиях, допускающих возможность оледенения. Запишем, следуя работе [1], одномерную нестационарную модель процессов в газе, транспортируемом по морскому газопроводу.

Уравнение неразрывности

Модель 1

1^ = * «>

уравнение движения

д(ри) д . , ulul

+ т-(р + ри-2) = -Xp-^ + pg cos a(z)-, (2)

dt dz ' ' ' 4R

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

d(pe) d ( ( p\ .

+ — I pu I e H— = — lo + pug cos a(z), (3)

dt ' dz\r V P

где

e = e + u2/2; (4)

калорическое уравнение

3 с

e = 1п(1 + ад; (5)

уравнение состояния Редлиха-Квонга

= hpT___ср2

р 1-Sp (1 + <V)TV2' [>

f f

h = Rg/M, M = Y, ni m, Y, ni = 1,

i=l i=l

S = b/M, с = a/M2, b = QbRg Tc/pc, a = Qa(Rg )2 T^'5/pc,

в котором h, с, S — размерные постоянные; Rg — универсальная газовая постоянная; mi,ni — молекулярный вес и доля г-й составляющей газовой смеси соответственно; f — количество компонент газовой смеси; cv — удельный коэффициент теплоемкости смеси газа того же состава, но в состоянии идеального газа; Qa, — числа, определяемые для заданного химического состава газовой смеси по значениям критических температуры Tc и давления pc согласно таблицам, приведенным в [2].

Система уравнений (1)-(6) дополняется начальными и граничными условиями, соответствующими рассматриваемой задаче.

В модели 1 независимыми термодинамическими переменными являются плотность р и температура T газа, через которые выражаются внутренняя (е) и полная (е) энергии (4). Поэтому в качестве уравнения состояния (6) в модели 1 используется уравнение Редлиха-Квонга — одно из наиболее точных двухпараметрических уравнений состояния в широком диапазоне изменений р, p и T вплоть до сверхвысоких давлений [2]. В работе [3] для уравнения состояния Редлиха-Квонга (6) получено калорическое уравнение (5). Для замыкания модели 1 необходимо найти выражение для слагаемого ш, входящего в уравнение энергии (3).

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

Величина ш выражается через qw — радиальную составляющую вектора плотности потока внутренней энергии (вектора потока тепла) на внутренней поверхности газопровода в z-м сечении:

/ ш dv = — ф q ■ n ds,

Sn

ш dv = шпЁ2 Sz,

n

где П — область, ограниченная поперечными сечениями газопровода, проходящими через г и г + 5г, и боковой поверхностью газопровода между данными сечениями. Тепловые условия на внешней поверхности газопровода на масштабах 5г допустимо

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

УС € [г, г + 5г], Чш(С,Л) = Чш(г,Л) ^

^ ^ Ч ' п За = J ! qw « 2п К 5 г qw (г,Л) ^

Бп 0 2

2 2ч- (г Л)

—> и;(г,1)'кН Зг « — 2тт К Зг ди1(г,1), ш(г^) « —-

К

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

Для установившегося варианта модели 1 в книге [4] и статье [5] получены явные выражения для потока тепла = ч(К) при наличии льда на боковой поверхности газопровода и при его отсутствии.

В ряде нестационарных задач при отсутствии слоя льда допустимо и в нестационарной модели 1 использовать формулу для ш, по форме аналогичную найденному в книге [4] стационарному выражению. Расчеты по модели 1 в этом приближении представлены в статье [1]. В настоящей работе исследуется погрешность такой замены. Для этого рассматривается нестационарная математическая модель 2 теплообмена между потоком газа в г-м сечении и окружающей средой для газопровода, имеющего два слоя обшивки — внутренний, состоящий из стали, и внешний — из бетона. Положим, что изменения как в окружающей среде, так и в потоке газа вдоль г на масштабах 5г ~ К пренебрежимо малы (данная ситуация характерна для большинства практических задач). При этом все функции в модели 2 будут зависеть от г, Л и, кроме того, параметрически от г через зависимость от г температуры газа Т(г,Л). Положим, что оледенение отсутствует. Примем, что поведение температуры газа во времени в г-м сечении является известной функцией. Выберем характерное поведение температуры газа в сечении газопровода, соответствующее расчетам [1]. Оно имеет следующий вид:

т = Ъ = а = Ъ(Т* — Тос), с = Тоо,

где Ти — заданная температура газа в заданный безразмерный момент времени Ли; Тж — минимальная температура газа.

Модель 2

дТ

= А1ВД), ге(Д,Й1), *€(0,*л); (7)

Л = 0, Тг(г)= Т *; (8)

Л € (0,£), г = К, Тг(К,Л)= Т(г,Л); (9)

г е (о, г), дТ2

дТ\

Я\, Т\=Т2, Ах—— = Л2

Р2С2-

дг

г е (о,г),

ру СУ

дЬ

Х2Ь(Т2), г = о,

г = я2,

Ху ЦП),

г = о,

дт2

дг 2 дг

г е (Й1,Й2), г е (о,г); Т2(г)= Т *;

Т2 — Ту, Х2 —— — л.

дг

дТу

Г е (Й2,Й2 + 6*), Ту (г) = Т*;

дг

г е (0,1]

г е (о,г), г = Я2 + 6*, Ту = Т*.

11)

12)

13)

14)

15)

16)

Модель 2 легко обобщается на большее число слоев обшивки газопровода.

В ней уравнения (7), (11) и (14) — одномерные линейные уравнения теплопроводности в слоях стали, бетона и эффективном тепловом погранслое воды. В этой модели (14)—(16) — совокупное влияние сложных процессов вокруг внешней поверхности газопровода, таких как немонотонное поведение плотности морской воды в окрестности температуры фазового перехода Т* , вызванные этим конвективные течения, донные течения, а также неосесимметричность условий контакта газопровода с донным грунтом учитываются с помощью введения эффективного параметра 6*, имеющего смысл толщины эффективного теплового погранслоя на внешней поверхности газопровода, в пределах которого передача тепла описывается линейным уравнением теплопроводности (14). Информация о 6* может быть получена из решения обратной задачи, которая базируется на экспериментальных данных в условиях, близких к реальным. Оценка возможной величины 6* в установившихся режимах приведена в книге [4] и статье [5].

Выражения (8), (12), (15) — начальные данные, которые соответствуют условию теплового равновесия в начальный момент времени, модель 2 легко обобщается на другие начальные данные; (10), (13) — условия равенства температур и тепловых потоков на границах областей. В рамках модели 2 поведение температуры газа в граничном условии (9) считается известным:

Т (г, г) = Т (г).

В частности, Т(г) может быть задана законом, приведенным выше. Расчет по модели 2 позволяет найти ш(г, г):

}(г,г) = -

2д(К,г) _ 2 дТ1

Я

дг

(17)

(18)

Модель 2 должна входить составной частью в модель 1. В общем случае решение системы уравнений модели 1 не расщепляется, ш(г,г) в (18) зависит от Т(г, г), являясь функцией (г, Т(г,г), Я¿, 6^, Х^, р^, е^, 6*, Т*). Решение по общей модели 1 выходит за рамки настоящей статьи.

Модель 2 предполагает отсутствие слоя льда на внешней поверхности газопровода. Момент времени I, начиная с которого температура Т21 на внешней поверхности газопровода опускается до температуры Т*, определяется из расчетов по модели 2. Условие Т2|Й2 ^ Т* необходимо (но не достаточно) для начала процесса оледенения.

к

г

При наличии слоя льда толщиной 6((г) модель теплообмена видоизменяется (см. модель 3).

Расчет по модели 2 преследует следующие цели:

1) определение момента времени г начала оледенения, вплоть до которого модель 2 имеет смысл;

2) оценка отклонений нестационарных потоков тепла д(Я,г) и д(Я2,г) от их значений до(Я) и до(Я2), рассчитанных по стационарному варианту модели 2. Величины отклонений Ад(Я,г) = до (Я) — д(Я,г) и Ад(Я2 ,г) = до (Я2) — д(Я2,г) позволяют оценить погрешность замены нестационарной модели 2 ее стационарным вариантом при расчете значения источника (стока) ш в модели 1. В случае допустимости такой замены интегрирование уравнений модели 1 существенно упрощается, так как параметр ш явно выражается через Т(г, г) и другие параметры задачи.

Решение системы уравнений модели 2 в стационарном варианте. При условии ■§[ = 0 уравнения теплопроводности (7), (11) и (14) легко интегрируются:

Щг)= А + Е11п г, г = 1, 2,У. (19)

Из граничных условий (10), (13), (16), (17) определим константы Ai, Bi

T * - T0

Bi =

Ri Ai R2 Ai R2 + S* ln — + — ln ——b — ln ■

R A2 Ri Av R2

Ai = To - Bi ln R, B2 = (Aí/A2)Bí, Bv = (Ai/Av)B1, (20)

Av = T* - Bv ln(R2 + S*), A2 = Av + (Bv - B2) lnR2.

Здесь T0 — некоторая заданная температура газа. Выражение для потоков до(R) = и </0(^2) = А2^-|Д2 находим из решений (19), (20) в явной форме:

AiBi /и \ A2B2 ab ,01ч

qo(R) = -s-, qo(R2) = = ■ (21)

R R2 R2

Система уравнений нестационарной модели 2 решалась численно методом сеток. Вводились характерные значения rx, tx переменных r, t, безразмерные переменные r, t выражались по формулам Г = r/rx, t = t/tx. (Далее штрихи у безразмерных переменных опущены.)

Интегрирование системы уравнений модели 2 не представляет трудности и может быть проведено с использованием как явных, так и неявных разностных схем. Для варианта явной разностной схемы применялась монотонная схема Самарского [6], аппроксимирующая задачу с первым порядком точности по At и вторым по Дг. Аппроксимация пространственных производных в граничных узлах имела первый порядок точности по Дг. Уменьшение шага Дг в 1.5 раза по сравнению с его значением Дг = 0.1 практически не вызвало в приведенном далее примере расчета изменения вычисляемых характеристик, поэтому дальнейшие расчеты проводились с безразмерным шагом Дг = 0.1. Для одномерных задач в цилиндрической системе координат (r, у, z) величины безразмерных шагов по времени At и по радиусу Дг, обеспечивающие устойчивость явной схемы, как известно [6], должны удовлетворять неравенству

At 1

<

J'm

Аг2 ^ 4ат

в котором ат для модели 2 является максимальным безразмерным коэффициентом температуропроводности, равным

ат = тах(аь а2, ау),

АгЛх

ргсг ГХ

% = 1, 2,У.

Результаты расчета по модели 2. Для тестовой задачи зададим следующие параметры (все величины указаны в системе СИ):

К = 0.51, К1 = К + 51 =0.55, К2 = К1 + 5Х =0.67,

5* =0.019, А1 = 24, Л2 = 1.7, Ау =0.6, р1 = 104,

р2 = 2300, ру = 1005, с1 = 450, с2 = 924, су = 4200,

Т* = 271, Т* = 272, Тк = 270, Тж = 266. (22)

В соответствии с законом Т(Л) при выбранных значениях Тк, Тж, Лк температура газа за сутки опускается от первоначальной величины Т* до температуры 268.2 К. Низкая температура фазового перехода в (22), равная 271 К, характерна для концентрации соли в морской воде 37% [7], а низкая температура воды (Т* = 272 К) — для ряда районов, например Баренцева моря.

Из проведенных расчетов тестовой задачи следует, что время I, начиная с которого выполняется неравенство Т2(К2,Л) ^ Т*, составляет I = 19.14 ч. На рисунке представлены распределения температуры в обшивке газопровода и эффективном тепловом погранслое воды для моментов времени Л1 = 2 ч, Л2 = 10 ч, Лз = 19.4 ч.

Г, К

272.0

271.0 -

270.0 "

269.0 -

268,

66 70

г, см

Распределение температуры в обшивке газопровода и эффективном тепловом погранслое воды для разных моментов времени 1 — I = 2 ч; 2 — I =10 ч; 3 — I = 19.14 ч.

Приведем рассчитанные размерные значения тепловых потоков с шагом в 2 ч на внутренней = и внешней = Аг^2-^ поверхностях газо-

провода:

г, ч.......................

д(й,4),Дж/(с.м2) д(Я.2,г), Дж/(е-м2)

4 6 8

0.425 0.433 0.442 0.095 0.148 0.191

10 12 14 0.452 0.463 0.474 0.225 0.253 0.276

16 18 19.4 0.485 0.494 0.50 0.296 0.313 0.322

Установившиеся значения потоков до (Я) и до (Я2), рассчитанные по формулам (19)—(21) для параметров (22) при температуре газа Т0 = Т({) = 268.31, равны 0.436 и 0.332 соответственно.

Выводы из расчетов по модели 2. Вывод 1. В тестовой задаче величина отклонения Ад(Я, Ь) = до (Я) — д(Я, Ь) на внутренней поверхности газопровода в течение всего процесса теплообмена вплоть до момента I не превышает нескольких процентов от стационарного значения до (Я).

Вывод 2. Величина отклонения Ад(Я2,Ь) = до(Я2) — д(Я2, Ь) потоков на внешней поверхности газопровода уменьшается с течением времени и в момент I возможного начала оледенения составляет 3% от стационарного значения до(Я2).

Вывод 3. Стационарный поток до(Я2) на внешней поверхности газопровода в течение всего процесса больше, чем нестационарный поток до(Я2,Ь).

Качественно выводы 1-3 предсказуемы, однако для каждого набора параметров (22) задачи количественная оценка погрешностей Ад(Я,Ь) и Ад(Я2,Ь) может быть найдена только в результате расчета.

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

Модель 3

дТ

= А1ВД), ге(Я,Я 1), ¿>£; (23)

Ь = г, Т1(г) = то (г); (24)

ь е (0,г), г = Я, тг(Я,г) = т(г,Ь); (25)

дТ дТ

*>*> г = Дь П=Т2, (26)

дТ2

Р2С2—- = Х2Ь(Т2), ге(ЯиЯ2), £ > (27)

дЬ

г = I, Т2(г)= Т°(г); (28)

дТ дТ

*>*> Г = Д2, Т2=Ти (29)

дТ

Р1С1-± = \,МТ,Х ге(я2,я2 + б1(г)), ¿>£; (30)

дЬ

дТ

¿>£, я, = я2 + Ш, А

Ь = г ё^ ) = 0;

дТу

— А^; —-

п1 дг

= да»

t>t, Ri = R2 + Si(t), Ti(Ri ) = П; (32) dTv

PvCv~z\- = \L(TV), r e (Ri,Ri + S*i), t > t; (33) dt

t = t, Tv (r) = T0(r); (34)

t>t, r = Ri, Tv = (35)

t>i, r = Ri + S*i, Tv = T *. (36)

В модель 3 наряду с уравнениями теплопроводности (23), (27), (30), (33) в слоях стали, бетона, льда и эффективного теплового погранслоя воды входит условие Стефана (31), выражающее тепловой баланс на фронте фазового перехода вода-лед. Оно служит уравнением для определения изменяющейся толщины слоя льда Si (t). Наличие изменяющегося во времени слоя льда существенно усложняет решение уравнений модели 3.

Задача Стефана исследовалась многими авторами. Существует обширная литература, посвященная как вопросам существования и единственности решения задачи Стефана (например, [8]), так и расчету теоретических и эмпирических зависимостей толщины слоя льда Si от всевозможных факторов (в частности, [7, 9]). В модели (33)-(36) толщина Sti эффективного теплового погранслоя воды в общем случае может отличаться от соответствующей величины S* в модели 2, поскольку, как известно [7], перед фронтом кристаллизации формируется слой более соленой морской воды, влияющий и на значение T*, и на процессы вокруг оледеневающего газопровода. В модель 3 может быть включен учет этих факторов. Интегрирование системы уравнений модели 3 будет представлено в наших дальнейших публикациях. Здесь ограничимся исследованием допустимых вариантов упрощения модели 3. Как следует из расчетов по модели 2 (вывод 2), на момент t начала оледенения погрешность замены нестационарного потока q(R2,t) стационарным qo(R2) невелика. Это позволяет предположить, что и при наличии медленно изменяющегося слоя льда значения нестационарного q(Ri,t) и стационарного q(Ri) потоков будут близки. Для проверки этой гипотезы рассмотрим упрощенную модель теплообмена газа с окружающей средой при наличии слоя льда (модель 4).

Модель 4

В ней в отличие от модели 3 считается, что слой льда заданной толщины Sio в момент времени t = to сформирован и остается далее неизменным на интервале времени Dt. Оценка допустимости предположения Si(t) = Sio, t G (to,to + Dt), приведена ниже (вывод 5). Начальные распределения температур To(r), T°(r), To(r) при t = to в модели 4 задаются из решения соответствующей стационарной задачи при температуре газа T(t) = const = Toi, Toi < T(i) = 268.31.

В модель 4 не входят уравнения (33)-(36) модели 3, поскольку распределение температуры Tv (r, t) в модели 3 влияет на остальные характеристики только через величину Si (t), а она в модели 4 считается постоянной и известной. Температура газа To на всем интервале времени Dt также считается постоянной и меньшей, чем температура газа Toi, при которой рассчитываются стационарные распределения температур Ti°(r), T°(r), To(r). С учетом указанных изменений уравнения (23)-(29) переносятся в модель 4. Уравнения (30)-(32) заменяются следующими:

dT

Pici-^г = WTi), г е {R2,Ri = const = Д2 + <ы t G (to,to + Dt)] (37) dt

г = г0, Щт) = Т?(т); (38)

г е (¿0,го + Щ, т = Щ = К2 + Зю, Т1(Е1 ) = Т*. (39)

Для численного интегрирования системы уравнений (23)—(29), (37)—(39) модели 4 использовался тот же алгоритм, что и для численного интегрирования уравнений модели 2. Для задания начальных условий в модели 4 найдем ее решение в стационарном варианте.

Решение системы уравнений модели 4 в стационарном варианте. Как

и в стационарном варианте модели 2 ((19),(20)), распределения температур Т^(т) (г = 1, 2,1) в слоях стали, бетона и льда имеют вид (19), константы Л^, ВI определяются граничными условиями модели 4 и выражаются через параметры задачи по формулам

Т* — То1

Вл

Й1 А1 Я2 А1 Я2 + Зю 1п — + — 1п ——Ь — 1п ■

Е А2 Е1 А1 Е2

Л1 = Т01 — В11п Е, В2 = (А1/А2)В1, В1 = (А1 /А{)ВЪ (40)

Л = Т* — В11п(Е2 + Зю), Л 2 = Л1 + (В1 — В2)1пЕ1.

Результаты расчетов по модели 4. Для тестовой задачи модели 4 выберем те же значения параметров, при которых рассчитывалась тестовая задача модели 2. Величины остальных параметров зададим в системе СИ:

Т01 = 267, Т0 = 265, З10 = 0.03, р1 = 928,

С = 2100, А1 = 2.3, т = 10 800, д = 335 000. (41)

В тестовой задаче считается, что температура газа в течение всего интервала времени тг = 3 ч поддерживается постоянной, равной 265 К. Кроме того, в течение 3 ч постоянной является толщина Зю слоя льда 3 см. Приведем результаты расчета по модели 4 нестационарных потоков тепла д(Е1 ,г) со стороны газопровода на фронте фазового перехода (т = Е1):

г, ч.......................... 123

я(Щ,г), Дж/(с • м2) 0.508 0.591 0.613

Стационарное значение потока д(Е1), рассчитанное по формулам (19), (40) при температуре газа 265 К, равно

дЕ) = 0.619.

Проведенные расчеты по модели 4 приводят к следующему выводу.

Вывод 4• В исследуемой тестовой задаче стационарный тепловой поток А;^- |й к границе фазового перехода Е1 со стороны льда в любой момент времени из интервала больше, чем нестационарный поток А; 1 й (как и следовало ожидать из физических соображений). Расчет по модели 4 позволяет оценить отклонение д(Е1) — д(Е1 ,г). В рассмотренной задаче оно на протяжении интервала времени Вг было незначительным.

Оценка максимальной величины приращения АЗ1 слоя льда за время Вг. Рассчитаем величину приращения слоя льда АЗ 1 за время Вг

АЗ1 = З1(г0 + т) — з^ц)

из уравнения Стефана (31), упростив его следующим образом: заменим нестационарный поток Л; 1 на стационарный поток А; 1 , тем самым увеличив количество образовавшегося за ВЬ льда, и соответственно значение ДЗ;, поскольку > Ví € (вывод 4). Величина рассчитывается по при-

веденному стационарному решению (40) и выражается через ДЗ; так:

Аг —

dr

А'(Г* Тш) R, = R, + S„ + AS,. (42)

Увеличим еще больше ДЗ;, положив в уравнении Стефана (31) тепловой поток

дТ„ I аг |д,

^Т^Н и со стороны воды к границе фазового перехода равным нулю.

Примем, что скорость ^ на интервале времени to < t < to + Dt можно считать постоянной и соответственно функцию 3¡(t) — линейной. Это оправдано, если окажется, что A3¡ — S¡o на интервале Dt:

Si(t) = ôl0 + ^t. (43)

Правая часть уравнения Стефана (31) с учетом этих упрощений запишется в виде

уравнение (31) с учетом (42) и (44) преобразуется в следующее трансцендантное уравнение относительно A3¡ :

■п* I + + + (45)

R Х2 Ri X¡ R2

Решение уравнения (45) при значениях параметров (22), (41) дает величину скорости оледенения

^ = 0.0007 м/ч. dt

Она является гипотетической максимальной, реальное значение скорости приращения льда ^ в этих условиях меньше.

Вывод 5. Максимально возможное в рассматриваемой задаче приращение льда A3¡ составляет 7% от начального значения S¡o = 0.03 м. Это, с одной стороны, свидетельствует о допустимости выбора линейной зависимости (43) на этом интервале времени, с другой — служит обоснованием реалистичности модели 4.

Литература

1. Груничева Е. В., Курбатова Г. И., Попова Е. А. Нестационарное неизотермическое течение смеси газов по морским газопроводам // Матем. моделирование. 2011. Т. 23, № 4. C. 141—153.

2. Рид Р., Праусниц Дж., Шервуд Т. Свойства газов и жидкостей / пер. с англ.; под ред. Б. И. Соколова. 3-е изд., перераб. и доп. Л.: Химия, Ленингр. отд., 1982. 592 с. (Reid Robert C., Prausnitz John M., Sherwood Thomas K. The properties of gases and liquids.)

3. Ермолаева Н. Н., Курбатова Г. И. Анализ подходов к моделированию термодинамических процессов в газах при высоких давлениях // Вестн. С.-Петерб. ун-та. Сер. 10. Прикладная математика. Информатика. Процессы управления. 2013. Вып. 1. С. 35—45.

4. Курбатова Г. И., Попова Е. А., Филиппов Б. В. и др. Модели морских газопроводов. СПб.: С.-Петерб. гос. ун-т, 2005. 156 с.

5. Филиппов В. Б. Модель морского газопровода с учетом оледенения // Вестн. С.-Петерб. ун-та. Сер. 1. Математика. Механика. Астрономия. 2004. Вып. 1. С. 103—111.

6. Самарский А. А., Гулин А. В. Устойчивость разностных схем. М.: Едиториал УРСС, 2005. 384 с.

7. Доронин Ю. П., Хейсин Д. Е. Морской лёд. Л.: Гидрометеоиздат, 1975. 320 с.

8. Мейрманов А. М. Задача Стефана. Новосибирск: Наука, 1986. 240 с.

9. Лейбензон Л. С. Собрание трудов: в 4-х т. М.: Изд-во АН СССР, 1955. T. 3. 678 с.

References

1. Grunicheva E. V., Kurbatova G. I., Popova E. A. Nestatsionarnoe neizotermicheskoe techenie smesi gazov po morskim gazoprovodam [Nonstationary nonisothermal flow of gas mix in offshore gas pipelines]. Matem. modelirovanie [Mathematical models and Computer simulations], 2011, vol. 23, no. 4, pp. 141—153. (In Russian)

2. Reid Robert C., Prausnitz John M., Sherwood Thomas K. Svoistva gazov i zhidkostei [The properties of gases and liquids]. Leningrad, Chemistry Publ., 1982, 592 p. (In Russian)

3. Ermolaeva N. N., Kurbatova G. I. Analiz podkhodov k modelirovaniiu termodinamicheskikh protsessov v gazakh vysokikh davleniiakh [The analysis of the approaches to the modeling of thermodynamic processes in gas flow at hyperpressure]. Vestnik of St. Petersburg University. Series 10. Applied Mathematics. Computer science. Control processes, 2013, issue 1, pp. 35—45. (In Russian)

4. Kurbatova G. I., Popova E. A., Filippov B. V. e. a. Modeli morskikh gazoprovodov [Models of sea gas-pipelines]. St. Petersburg, St. Petersburg University Press, 2005, 156 p. (In Russian)

5. Filippov V. B. Model' morskogo gazoprovoda s uchetom oledeneniia [Sea gas pipe line model with surface freezing effect consideration]. Vestnik of St. Petersburg University. Series 1. Mathematics. Mechanics. Astronomy, 2004, issue 1, pp. 103—111. (In Russian)

6. Samarskii A. A., Gulin A. V. Ustoichivost' raznostnykh skhem [Stability of difference schemes]. Moscow, Editorial URSS Publ., 2005, 384 p. (In Russian)

7. Doronin Iu. P., Kheisin D. E. Morskoi led [Sea ice]. Leningrad, Gidrometeoizdat Publ., 1975, 320 p. (In Russian)

8. Meirmanov A. M. Zadacha Stefana [The Stefan problem]. Novosibirsk, Nauka Publ., 1986, 240 p. (In Russian)

9. Leibenzon L. S. Sobranie trudov: v 4-kh t. [Collected papers: in 4 vol.]. Moscow, AS USSR Publ., 1955, vol. 3,678 p. (In Russian)

Статья рекомендована к печати проф. Н. В. Егоровым. Статья поступила в редакцию 30 апреля 2015 г.