CC BY
47
Вычислительные технологии
Том 11, № 6, 2006
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДВУМЕРНЫХ
НЕСТАЦИОНАРНЫХ ТЕЧЕНИЙ ГАЗА ЧЕРЕЗ ПОРИСТЫЕ ТЕПЛОВЫДЕЛЯЮЩИЕ ЭЛЕМЕНТЫ*
В. А. Левин, Н.А. Луценко Институт автоматики и процессов управления ДВО РАН,
Владивосток, Россия e-mail: levin@iacp.dvo.ru, nickl@inbox.ru
The numerical method, based on the combination of explicit and implicit finite difference schemes, is proposed for modeling the time-dependent two-dimensional gas flow through porous heat-evolutional elements, which may arise from natural or human related disasters (as destroyed block of Chernobyl atomic power station). The cooling of convergent porous elements and step changed porous elements is investigated. A cooling process is shown to be affected by the shape of the heat-evolutional element.
Введение
В настоящей работе рассматривается течение газа через твердый пористый тепловыделяющий элемент. Такое движение может возникать, например, при охлаждении очагов тепловыделения в пористых средах, которые могут появляться в результате природных или техногенных катастроф. Подобная модель предложена для описания процесса охлаждения разрушенного блока Чернобыльской АЭС [1]. Исследование рассматриваемых течений газа может быть полезным при моделировании новых технологических процессов (например, в металлургии и химической промышленности), создании новых охлаждающих устройств в электронной промышленности.
Стационарный режим воздушного охлаждения при естественной конвекции пористых тепловыделяющих элементов исследован в [1]. Отличительной особенностью модели является открытость саморазогревающейся пористой среды в атмосферу снизу и сверху. И хотя применяемые уравнения являются классическими и использованы в той или иной модификации во многих работах по теории фильтрации, новый тип краевой задачи для них, возникший при анализе конкретных условий охлаждения аварийного энергоблока Чернобыльской АЭС, привел к открытию новых физических эффектов. В [2, 3] рассмотрены стационарные режимы охлаждения пористых тепловыделяющих элементов принудительно нагнетаемым потоком газа. В [4, 5] предложен численный метод для моделирования нестационарных одномерных течений газа через пористые элементы с выделением тепла.
* Работа выполнена при финансовой поддержке Президента РФ (грант № МК-6493.2006.1), Российского фонда фундаментальных исследований (грант № 06-01-96020-р_восток_а), ДВО РАН (проект № 06-III-В-03-079).
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.
В [5, 6] показано, что при моделировании движения газа через пористую тепловыделяющую среду необходимо учитывать температурную зависимость его вязкости. В [5] обнаружено явление неограниченного разогрева и расплавления пористого тепловыделяющего элемента при докритических краевых условиях, при которых возможно существование установившегося стационарного режима.
Настоящая работа посвящена исследованию плоских (двумерных) нестационарных течений газа через пористые тепловыделяющие элементы. Описывается численный метод, разработанный для моделирования таких процессов. Решается задача о включении принудительной фильтрации в момент начала тепловыделения в твердой фазе для тепловыделяющих элементов плавно сужающейся и ступенчато сужающейся формы.
1. Математическая модель
Пористый тепловыделяющий элемент высотой Н предполагается однородным и неподвижным, с боков он ограничен нетеплопроводными стенками, а сверху и снизу открыт. В его нижнюю часть под давлением подается холодный газ, который движется снизу вверх через пористую среду, нагреваясь в результате теплообмена, и вытекает в свободное пространство с заданным давлением. Модель строится в рамках модели двух взаимодействующих континуумов [7]. Для описания динамики газа будем использовать уравнение сохранения импульса для пористых сред, которое является более корректным, чем классическое уравнение Дарси и может применяться при достаточно больших числах Рейнольдса. Будем считать, что процесс тепловыделения в твердой фазе происходит в результате протекания "химической" реакции, тепловыделение прямо пропорционально концентрации реагирующего вещества С, скорость убывания которой прямо пропорциональна самой концентрации, объем и масса конденсированной фазы изменяются незначительно, и этими изменениями можно пренебречь. Предположим, что интенсивность межфазового теплообмена пропорциональна разности фазовых температур в рассматриваемой точке среды, для газа справедливо уравнение состояния совершенного газа. Учитывая в уравнениях энергии фаз теплопроводность и работу внутренних сил в газе, запишем систему уравнений, моделирующую нестационарное течение газа через пористый тепловыделяющий элемент:
дТ
(1 - а) РсСс= -а (Т - Т9) + Оо (1 - а) С + (1 - а) А ^
д 2Т дх2'
арср
дТ
дг
+ дТЛ Т (д'р др\ +2 + д / дТА
+^ зх,) =а (Т - Т)+аи+% дхг) + ак~Х +адж О дхг)
Р (1 + X (1 - а))(+ ^-
дг
дх,
Р = рЯТд
др дх
РЯг - а—Удг,
f = -к2С
V
дар 1 даргд ^
дг
+
дх,-
(1)
Здесь и далее а — пористость; сс — теплоемкость конденсированной фазы; ср — теплоемкость газа при постоянном давлении; д — ускорение силы тяжести; кг — коэффициент проницаемости конденсированной фазы; к2 — коэффициент, определяющий уменьшение тепловыделения; р — давление газа; О0 — константа, определяющая интенсивность тепловыделения; Я — газовая постоянная; Ь — время; Т — температура конденсированной
фазы; Тд — температура газа; х — эйлерова координата; уд — скорость газа; а — константа, определяющая интенсивность межфазового теплообмена; Л — теплопроводность конденсированной фазы; Лд — теплопроводность газа; / — динамическая вязкость газа; р — плотность газа; рс — плотность конденсированной фазы; х — коэффициент присоединенной массы, учитывающий инерционное взаимодействие фаз при их ускоренном относительном движении. Индексы г, ] — номера декартовых координат (для двумерного случая: 1 — горизонтальная, 2 — вертикальная), по повторяющимся индексам производится суммирование.
При исследовании неизотермической фильтрации жидкости динамическая вязкость, как правило, полагается зависящей от температуры [8], однако при моделировании неизотермической фильтрации газа вязкость часто принимается постоянной. В [5, 6] показано, что при моделировании движения газа через пористую тепловыделяющую среду необходимо учитывать температурную зависимость вязкости газа. Далее будем считать, что динамическая вязкость газа зависит от температуры по формуле Сазерленда, а теплопроводность газа пропорциональна его вязкости:
Т1-5 с
// — С д л — Ср //
^ — Сз 1--77т , Лд — Б"
Св2 + Тд Рг
где Рг — число Прандтля для газа. Введем обозначение: и — аьд — скорость фильтрации газа. Из последнего уравнения системы (1) найдем выражение для концентрации реагирующего вещества: С — ехр(-к2£). Далее введем безразмерные переменные следующим образом: х — Нх, £ — и — и*й, где £* и и* — характерные значения времени и скорости фильтрации газа; р — р*р, р — р*р, Т — Т*Т, Тд — Т*Тд, где р*, р*, Т* — давление, плотность, температура газа при "нормальных" условиях. Будем использовать следующие параметры подобия:
8Ь — инН*, 811 — -аН-) 812 — Ей — рри_, Ее —
Н —СССи* р*Сри* р*и* СрТ *
Ке — риН, Рв1 — и*-сСсН, Рв2 — Рг, Рг — ,
Сз1л/Т* Л Сз1л/Т* дН
к1 Я0£* , , ~ Сз2
п — Н2, Я — -7Г, е — 2^ Сз2 — Т^-
Н рсСсТ* Т*
Заметим, что так как теплопроводность газа мала, Ре2 >> 1 и слагаемым, его содержащим, можно пренебречь. Для случая, когда все параметры состояния модели зависят от двух декартовых координат, перепишем систему (1) в безразмерных переменных, опустив тильду:
дТ 8Ь811 Т Я ( д2Т
(Т - Тд) + Я ехР (-е£) + РГХ^
д£ 1 - а ^ д' ' ^ к ^ У ' Ре^ дх2 '
г '
дТд 8Ь дТд N 8Ь812 ^ ^ ^ /др 8Ь др N 8ЬЕе ^ 2 Т1-5
\ Л I I --I / I Л I II \ Л ^
р -ттт + — и^ —-2 (Т - Тд) + ЕиЕе —у + — иг^~ + -— > и
д£ а г дхг) а д V д£ а г дхг У аИ,еп г сз2 + Тд'
1 + х (1 - а) + и.диЛ — -Еи8ь- и1 Т1-5
а р \ д£ а г дхг / дх1 И,еп 1 сз2 + Тд,
1 + X (1 - а) (ди 8Ь ди2 N С1 др 8Ь 8Ь ТЙ1-5 -р т;--1--иг—— — -Еи8п---— р - —— и2:
д£ а дхг / дх2 Рг И,еп сз2 + Тд
dp Sh dpui о dt a dxi '
P = PTg
На входе в пористый элемент известны температура газа и давление. На выходе известно давление, так как истечение газа происходит в открытое пространство. Известны также условия теплообмена на входе и выходе из пористого элемента и на ограничивающих непроницаемых стенках. Отличительной особенностью модели является то, что расход и скорость фильтрации газа на входе в пористый элемент неизвестны и должны определяться при решении задачи. Таким образом, краевые условия для системы (2) имеют следующий вид:
PL=o = Po (t) , Tg1
g |Х2=0
Tg0 (t)
dT
dx2
X2=0
Bi (T 1x2=0 - Tgo) , U1|X2=0
0,
P|X2 = 1 = Ph,
0Tn
dx2
X2 = 1
dT dn
xeG
dT
dT
0, g
dn
dx2
xeG
Bi (Tg 1x2=1 - TIx2=i)
x2=1
= 0, UnlxeG = 0-
(3)
Здесь Б1 = /ЗИ/Х, где в — коэффициент теплоотдачи; С — поверхность боковых стен; п — нормаль к поверхности боковых стен.
Для решения системы (2) необходимо также задать значения искомых величин в начальный момент времени.
0
2. Численный метод
Система уравнений (2), моделирующая нестационарное движение газа через пористый тепловыделяющий элемент, является нелинейной смешанной гиперболически-параболической системой уравнений и в общем случае не может быть решена аналитически. Для ее решения разработан численный метод, основанный на комбинации явных и неявных конечно-разностных схем. Уравнения энергии и сохранения импульса преобразуются в явные конечно-разностные уравнения, из которых определяются температура газа, температура твердой среды и скорость фильтрации газа. Уравнение неразрывности преобразуется в неявное конечно-разностное уравнение, из которого с учетом уравнения состояния совершенного газа определяется давление газа. Плотность газа находится из уравнения состояния совершенного газа.
Рассмотрим равномерную сетку с шагом h по пространству и шагом по времени т, равным rh2, где r = const. Нижние индексы при искомых сеточных функциях будут обозначать продвижение по пространству: m — по горизонтали, l — по вертикали; верхний индекс n — продвижение по времени. Тогда можно записать следующую систему конечно-разностных уравнений, аппроксимирующую исходную систему (2):
Tn+1 = Tn (1 _ ShSt1T ) i Sh r (rpn 4Tn + Tn + Tn + Tn \ +
Tm,l = Tm,l I 1 i _ a J + pe1 ' \1m+1,l 41m,l + 1 m-1,l + 1 m,l+1 + 1m,l-1) +
ShSt1r n
+ 1 _ Tgm,l + TQ exP (-£nT)
тп +1 _ Тп
1 дт,1 1 дт,1
8Ь812Л ЭМ (, 3Тп 4Тп +Тп А_
1 п I 2а и1т,1 \±01дт,1 + 41дт^1,1 ± Тдт^2,1)
аРт,1
ЭЪтН ~2а
-и.
2т, I
(±3Тдтт,г Т 4Тдт,гТ1 ± Т'Пт,1^ +
арт,1
Тп +
Тт,1 +
+
ЕиЕс
п
Рт,1
Рп РтЛ
Рп-1 +
Рт,I +
~2а
и
1т,1 (Р т+1,1 Рт-1,1) +
~2а
и
2т,I
Рт,1— 1 /
+
Тд
дт,^1,5 ЕсБЬг
5*2 + ТдПт,1
(и1т,г) + (и2т,г)
ип+1 _ ип I 1 _ и1т,1 _ и1т,1 I 1
и
2а ^1т+1,г М1т-1'^ (1 + х (1- а)) Иер^г сЙ + Тд аЕиЭЬтЛ,
аБЬт
(Т
п
дт,1
1,5
ЭЪтН ~2а
и2т,1 (и1т,1+1 и1т,1-1
п дт,1
2(1 + х (1-а)) р;
п
т,1
(Рт+1,г - ^т-1,0 - и3 - и4, (4)
ип+1 _ ип (1 (ип
и2т,1 _ и2т,1 ' 1 - I и
2а
2т,1+1 и2т,1-1
аБЬт
Тд
п
дт,1
1,5
(1 + х (1 - а)) Иерт.г ^2 + Т
п 1 дт,1
п п п
и и и
аЕиЭЬтЛ,
2а -1тМ-2т+1>1 "2т-1,и 2(1 + х (1 - а)) рпт1^т
Рт,1— 1 /
аБЬт
(1 + х (1 - а)) Ег
^5 -
Р2+1
БЪтк
(ЭЬтЛ,
2пТ п+1 2а1дт,И
рп+1 + 1т+т+рт+1
и
п+1 2т,1+1
Пп гт,1
БЪтк
рп+1
рт,1
(г! „,п _ Пп п \ п+1 _ 1 тЛ
\Рт+1,1а1т+1,1 Рт-1,1 а1т-1,1) и1 и2, рт,1 _ т™+1 '
дт,1
Система конечно-разностных уравнений (4) аппроксимирует исходную систему (2) с краевыми условиями (3) со вторым порядком точности по пространству и первым порядком по времени. Для схем четного порядка точности типично преобладание дисперсионной ошибки [9], поэтому для ее уменьшения в систему (4) введены демпфирующие члены и1,.'' , и6. Эти члены имеют четвертый порядок и не изменяют формальную точность метода:
и1 _ — (Рт+2,г - 4рт+1,1 + 6Рт,1 - 4РПт-1,1 + Рт-2,0 ' и* _ СОПй1'
Ш-,
и _ —1 (Рт,1+2- 4Рт,г+1+бРт,1- 4Рт,г-+рт,—), и* _ с°п^
^3
^2
ип — 4ип и 4и
,* V 1т+2,1 Шз
1т+1,1
+ Кт,! - 4и1т-1,1 + и1т-2,0 ' _ СОП^
и4 — з (и1т,г+2 4и1т,1+1 + 6и1т,1 4и1т,1-1 + и1т,1-2)
и,
и4* _ сопэ^
Ш5
и,
ил — Аип
* \а2т+2,1 4и2т+1,1
+ 6и2т,1 4и2т-1,1 + и2т-2,\)
и* _ с°пэ1,
+
1
= -1 Кт;1+2 - 4<п,1+1 + 6< г - 4иПт1-1 + < , Ч* = сопэ!.
Во втором уравнении системы (4) знак ± обусловлен возможным изменением направления движения газа: во второй скобке верхнее значение знака выбирается при положительной горизонтальной скорости фильтрации газа в рассматриваемой точке, а нижнее значение — при отрицательной; в третьей скобке верхнее значение знака выбирается при положительной вертикальной скорости фильтрации газа в рассматриваемой точке, а нижнее значение — при отрицательной.
Краевые условия для сеточных функций тривиально получаются из (3), но следует добавить искусственные условия, не уменьшающие порядок точности разностной схемы, для скоростей фильтрации газа на входе и выходе из элемента и для давления газа на поверхности боковых стенок. Необходимо также добавить фиктивные точки, расположенные вне физической границы, и по аналогии задать в них условия для давления, скоростей фильтрации и температуры газа. Следует сделать оговорку, что уравнения (3) справедливы для положительной вертикальной скорости фильтрации газа на входе в тепловыделяющий элемент (т.е. при и2|х2=0 > 0). При нарушении этого условия следует условие известной температуры газа на входе в элемент заменить на условие известного градиента температуры газа на входе в элемент: дТа/дх21х =0 = 0.
Алгоритм нахождения искомых величин на каждом временном слое следующий: фиксируем индекс "горизонтали" т; последовательно меняя индекс "вертикали" п, решаем первое — четвертое уравнения системы (4); затем методом прогонки решаем пятое уравнение и далее тривиально решаем шестое уравнение системы (4); после этого меняем индекс "горизонтали" т и повторяем процесс. Таким образом, задавая начальные условия и последовательно продвигаясь по временным слоям, можно определить искомые величины в требуемый момент времени.
Доказать аналитически сходимость системы (4) к (2) не представляется возможным. Следует заметить, что даже если бы удалось аналитически доказать устойчивость системы (4), то это не приводило бы к доказательству сходимости, так как система уравнений (2) является нелинейной. Экспериментальное исследование сходимости системы конечно-разностных уравнений (4) к системе (2), проведенное на ряде модельных задач, показало, что сходимость имеет место при некотором ограничении на г, причем с уменьшением шага Н допустимое значение г увеличивается. Полученное экспериментально ограничение на шаг по времени оказывается более жестким, чем ограничение, полученное аналитически из анализа устойчивости уравнений системы (4), проведенного по спектральному методу Неймана с учетом принципа замороженных коэффициентов [10].
Численный эксперимент показал, что точность предлагаемого метода зависит от геометрии пористого элемента. Если конфигурация элемента такова, что задача сводится к одномерной и при этом устанавливается стационарный режим, то решение можно получить численно-аналитическим способом [3] и сравнить его с результатами расчетов по предлагаемому методу. Такое сравнение показало практически полное совпадение решений даже на достаточно грубой сетке (при шаге к = 0.05 погрешность составляла менее одного процента). Для элементов других конфигураций исследование точности проводилось на последовательности сгущающихся сеток. Анализ показал, что для сетки с шагом к = 0.025 погрешность в основном не превышает нескольких процентов для элементов различной геометрии. Но исключение составляет точность нахождения локальных экстремумов скорости фильтрации, которые могут возникать, например, при огибании газом горизонтальных препятствий. Погрешность вычисления таких пиковых значений скоро-
сти на сетке с шагом h = 0.025 может составлять несколько десятков процентов. Это естественно — достаточно грубая сетка не в состоянии точно уловить экстремум при быстром изменении искомой величины. Однако если первоочередной интерес вызывает нахождение температур, давления или плотности, то для элементов таких форм, при которых существуют пиковые значения скорости, также можно использовать достаточно грубые сетки.
Система конечно-разностных уравнений (4) может применяться при моделировании разнообразных нестационарных плоских задач, возникающих при изучении течения газа через однородный пористый тепловыделяющий элемент. Достоинством предложенного численного метода является отсутствие необходимости решения сложных систем алгебраических уравнений, благодаря чему c использованием геометрической декомпозиции возможно легкое и эффективное распараллеливание алгоритма.
3. Охлаждение пористых элементов плавно сужающейся формы
Рассмотрим тепловыделяющий элемент, у которого две боковые стенки элемента вертикальны и параллельны, а другие две — им ортогональны, но имеют более сложную форму: в нижней части вертикальны, а в верхней части плавно сужаются (рис. 1). Можно считать, что в этом случае движение газа через пористый тепловыделяющий элемент является плоским и описывается системой уравнений (2) с краевыми условиями (3). Для простоты будем рассматривать элементы, у которых сужающиеся стенки имеют угол п/4 к горизонтальной плоскости. В этом случае при расчете можно использовать равномерную сетку и искомые значения на границах будут определяться тривиально.
Рассмотрим следующую задачу. Тепловыделение в твердой фазе до начального момента времени отсутствует, давление на входе в элемент и на выходе из него соответствует атмосферному давлению на заданных высотах, следовательно, движение воздуха в элементе отсутствует. В начальный момент времени начинается тепловыделение в твердой фазе и одновременно происходит быстрый рост давления газа на входе в элемент.
Далее, если не оговорено особо, будем рассматривать систему (2) со следующими безразмерными параметрами:
Sh = 0.1, Sti = 4.94 ■ 10-3, St2 = 8.33, Eu = 8.33 ■ 104, Ec = 3.33 ■ 10-6, Re = 4.75 ■ 105, Pe1 = 1.687 ■ 107, Fr =1.02 ■ 10-2, Bi = 83.33, п = 10-10, (5)
-i-►
u
Рис. 1. Плавно сужающийся пористый элемент.
д = 1.647 ■ 10-4, е = 10-7, с,2 = 0.368, а = 0.3, х = 0.5. Эти данные соответствуют следующим размерным величинам:
Н =10 м, ^ = 1 с, и* = 1 м/с, Т = 300 К, р* = 105 Па, р* = 1.2 кг/м3
Дж 1П3 Дж щ3 Дж
рс = 2.2 ■ 103 кг/м3, сс = 9.2 ■ 102^—, а = 103 ^—, ср =10
кг • К м3 • К • с кг • К
-6 кг „ , _1П-8„2 ^ _ 1 п5 Дж
св1 = 1.458 ■ 10-6-==, с^2 = 110.4К, к1 = 10-8 м2, д = 10 ч
3
мс
' ' ' м3 ■ с'
Дж 2 7 1 Дж
Л =1.2^-—, д = 9.8 м/с2, к2 = 10-7 -, в = 10-^—. м ■ К ■ с с м2 ■ К ■ с
Пусть давление на входе в элемент быстро возрастет до 1.5 и после достижения этого значения остается постоянным. Тогда краевые условия получаются из (3) при следующих значениях:
ро = 1.5, Тдо = 1, рн = 1. (6)
Вначале рассмотрим тепловыделяющий элемент, ширина которого в верхней части Ьн в два раза меньше, чем ширина Ь0 в нижней. Пусть Ь0 = 2, плавное сужение элемента начинается с высоты х2 = 0.5, наклонные стенки, таким образом, имеют угол п/4 к горизонтальной плоскости. В этом случае происходит неограниченный разогрев тепловыделяющего элемента, который неизбежно закончится расплавлением твердой фазы и нарушением описанного процесса охлаждения. На рис. 2 приведены искомые величины через 10 ч после начала процесса. Для лучшей визуализации графики давления и плотности газа развернуты на 180°. Как видно из рисунка, в средней части элемента проходящий воздух достаточно эффективно охлаждает конденсированную фазу, так как движется здесь более быстро. А возле боковых стенок, как наклонных, так и вертикальных, охлаждение не столь эффективное и наблюдается более сильный разогрев твердой среды — газ здесь более разреженный и движется вверх медленнее. Встречая на своем пути наклонные стенки, газ стремится их обогнуть, поэтому в таких местах наблюдаются повышенные давление и горизонтальная скорость фильтрации. Нагреваясь в результате движения по элементу, газ начинает хуже охлаждать твердую среду, поэтому температуры газа и конденсированной фазы возрастают по направлению к выходу из элемента, плотность и давление газа при этом убывают, а скорость фильтрации также увеличивается. Наибольший разогрев твердой среды происходит вблизи боковых стенок возле выхода из элемента. Также в этом месте наблюдаются пиковые значения скорости фильтрации газа: как вертикальной, так и горизонтальной. Потоки газа, огибающие наклонные стенки, больше не встречают сопротивления, поэтому резко устремляются из элемента.
Далее рассмотрим тепловыделяющий элемент, ширина которого в верхней части Ьн в четыре раза меньше, чем ширина Ь0 в нижней. Пусть Ь0 = 2, плавное сужение элемента начинается с высоты х2 = 0.25, наклонные стенки, таким образом, имеют угол п/4 к горизонтальной плоскости. В этом случае также происходит неограниченный разогрев тепловыделяющего элемента, но более сильное сужение приводит к ухудшению охлаждения: температура твердой фазы и газа возрастает как в средней части элемента, так и у боковых стенок. График давления более выпуклый — оно убывает медленнее вблизи входа в элемент и быстрее у выхода из элемента. А график плотности, наоборот, более вогнутый — вблизи входа в элемент убывает быстрее, а у выхода медленнее. Вертикальная скорость фильтрации вблизи входа в элемент незначительно уменьшается. Объяснение
Рис. 2. Распределение температуры твердой среды, температуры газа, давления, плотности газа, вертикальной и горизонтальной скорости фильтрации газа по плавно сужающемуся пористому элементу при Ь^ = 0.5Ьо.
обнаруженному следующее: основная часть газа проходит через среднюю часть элемента, где нет препятствий, но так как эта часть стала уже, уменьшился и суммарный расход газа, что привело к ухудшению охлаждения, росту температуры и изменению других искомых величин.
Следует заметить, что в случае отсутствия сужения тепловыделяющего элемента, т. е. когда все боковые стенки вертикальны, при параметрах (5) и краевых условиях (6) имеет место установление стационарного режима охлаждения (точнее, установление квазистационарного режима, так как интенсивность тепловыделения со временем медленно убывает из-за "выгорания" реагирущего вещества, что приводит к медленному изменению всех остальных параметров). На рис. 3, а показано изменение температуры твердой среды со временем в самых горячих точках как плавно сужающихся элементов, так и не сужающегося. Таким образом, видно, что сужение элемента в верхней части приводит к существенному ухудшению охлаждения, вплоть до нарушения стационарного режима и перехода к катастрофическому разогреву.
Однако и для плавно сужающихся элементов возможно установление стационарных режимов охлаждения. На рис. 3, б показано изменение температуры твердой среды со временем в самых горячих точках плавно сужающихся элементов, рассмотренных выше, но при интенсивности тепловыделения в 10 раз меньше. Остальные параметры подобия определяются (5), а краевые условия — (6). Как видно из рисунка, в этом случае процесс охлаждения стабилизируется. Графики распределения искомых величин по элементам имеют вид, значительно схожий с рис. 2, хотя значения самих величин существенно отличаются от рассмотренного выше случая: температура твердой фазы и газа значительно ниже; график давления менее выпуклый, а плотности — менее вогнутый; скорость фильтрации, как вертикальная, так и горизонтальная, выше.
Таким образом, можно сделать вывод, что при охлаждении плавно сужающихся пористых тепловыделяющих элементов возможны как установление стационарного режима, так и неограниченный разогрев, ведущий к расплавлению твердой фазы. Плавное сужение элемента в верхней части приводит к существенному ухудшению охлаждения: основная часть газа проходит через среднюю часть элемента, где нет препятствий, поэтому
Рис. 3. Изменение температуры твердой среды со временем в самых горячих точках плавно сужающихся пористых элементов при Q = 1.647 ■ 10-4 (а) и Q = 1.647 ■ 10-5 (б).
ее сужение приводит к уменьшению суммарного расхода газа и, соответственно, к росту температуры и изменению других искомых величин. Наибольший разогрев твердой среды происходит вблизи боковых стен возле выхода из элемента.
4. Охлаждение пористых элементов ступенчато сужающейся формы
Рассмотрим тепловыделяющий элемент, у которого две боковые стенки элемента вертикальны и параллельны, а другие две — им ортогональны, но имеют более сложную форму: нижняя и верхняя части вертикальны, но различаются по ширине, нижняя часть более широкая, чем верхняя; эти части соединены между собой горизонтальной непроницаемой стенкой (рис. 4). Можно считать, что в этом случае движение газа через пористый тепловыделяющий элемент является плоским и описывается системой уравнений (2) с краевыми условиями (3).
Рассмотрим задачу, аналогичную исследованной выше, для плавно сужающихся элементов. Тепловыделение в твердой фазе до начального момента времени отсутствует, давление на входе в элемент и выходе из него соответствует атмосферному давлению на заданных высотах, следовательно, движение воздуха в элементе отсутствует. В начальный момент времени начинается тепловыделение в твердой фазе и одновременно происходит быстрый рост давления газа на входе в элемент. Пусть параметры подобия определяются (5), а краевые условия — (6).
Вначале рассмотрим тепловыделяющий элемент, ширина которого в верхней части Lh в два раза меньше, чем ширина L0 в нижней. Пусть L0 = 1, сужение элемента имеет место на высоте x2 = 0.5. В этом случае происходит неограниченный разогрев тепловыделяющего элемента, который неизбежно закончится расплавлением твердой фазы и нарушением описанного процесса охлаждения. На рис. 5 приведены искомые величины через 10 ч после начала процесса. Для лучшей визуализации графики давления и плотности газа развернуты на 180°. Как видно из рисунка, в средней части элемента проходящий воздух достаточно эффективно охлаждает конденсированную фазу, так как движется здесь более быстро. А возле боковых стенок охлаждение не столь эффективное, и наблюдается более сильный разогрев твердой среды (газ здесь более разреженный и движется вверх медленнее). Нагреваясь в результате движения по элементу, газ начинает хуже охлаждать твердую среду, поэтому температуры газа и конденсированной фазы возрастают по направлению к выходу из элемента, плотность и давление газа при этом убывают. В средней части эле-
Рис. 4. Ступенчато сужающийся пористый элемент.
Рис. 5. Распределение температуры твердой среды, температуры газа, давления, плотности газа, вертикальной и горизонтальной скорости фильтрации газа по ступенчато сужающемуся пористому элементу при Ь^ = 0.5Ьо.
мента вертикальная скорость фильтрации также возрастает по направлению к выходу из элемента, в отличие от краев, где она сначала увеличивается, а потом падает до нуля — горизонтальные стенки не пропускают газ. Встречая на своем пути горизонтальные непроницаемые стенки, потоки газа стремятся обогнуть их, максимально усиливаясь у углов, образованных соединением этих стенок с верхними вертикальными стенками. В этом месте наблюдается пиковое значение горизонтальной скорости фильтрации газа. После того как потоки газа огибают эти углы, вертикальная скорость фильтрации имеет локальный пик. Наиболее слабое охлаждение — в углах, образованных соединением нижних вертикальных стенок с горизонтальными непроницаемыми стенками: здесь температура растет наиболее интенсивно и наблюдается максимальный разогрев твердой среды. Слабый отвод тепла из этих углов объясняется очень малыми скоростями фильтрации газа. Также в этих местах давление повышено по сравнению со средней частью элемента.
Далее рассмотрим тепловыделяющий элемент, ширина которого в верхней части Ь^ в четыре раза меньше, чем ширина Ь0 в нижней. Также Ь0 = 1, а элемент сужен на высоте х2 = 0.5. И в этом случае происходит неограниченный разогрев тепловыделяющего элемента, но, как и для плавно сужающегося элемента, более сильное сужение приводит к ухудшению охлаждения: температура твердой фазы и газа возрастает как в средней части элемента, так и у боковых стенок. График давления немного более выпуклый — оно убывает медленнее вблизи входа в элемент. А график плотности, наоборот, более вогнутый — вблизи входа в элемент убывает быстрее, а у выхода медленнее. Вертикальная скорость фильтрации уменьшается по всему элементу. Аналогично задаче о плавно сужающемся элементе обнаруженному можно дать следующее объяснение: основная часть газа проходит через среднюю часть элемента, где нет препятствий, но так как эта часть стала уже, уменьшился и суммарный расход газа, что привело к ухудшению охлаждения, росту температуры и изменению других искомых величин.
Как отмечалось выше, в случае отсутствия сужения тепловыделяющего элемента, т. е. когда все боковые стенки вертикальны, при параметрах (5) и краевых условиях (6) имеет место установление квазистационарного режима охлаждения. На рис. 6, а показано изменение температуры твердой среды со временем в самых горячих точках как ступенчато
Рис. 6. Изменение температуры твердой среды со временем в самых горячих точках ступенчато сужающихся пористых элементов при Q = 1.647 ■ 10-4 (а) и Q = 1.647 ■ 10-5 (б).
сужающихся элементов, так и не сужающегося. Таким образом, так же, как и при плавном сужении элемента, ступенчатое сужение элемента приводит к существенному ухудшению охлаждения, вплоть до нарушения стационарного режима и перехода к катастрофическому разогреву.
Однако и для ступенчато сужающихся элементов также возможно установление стационарных режимов охлаждения. На рис. 6, б показано изменение температуры твердой среды со временем в самых горячих точках ступенчато сужающихся элементов, рассмотренных выше, но при интенсивности тепловыделения в 10 раз меньше. Остальные параметры подобия соответствуют (5), а краевые условия — (6). Как видно из рисунка, в этом случае процесс охлаждения стабилизируется, как и для плавно сужающихся элементов. Графики распределения искомых величин по элементам имеют вид, значительно схожий с графиками, приведенными на рис. 5, хотя значения самих величин существенно отличаются от рассмотренного выше случая: температура твердой фазы и газа значительно ниже; график давления менее выпуклый, а плотности — менее вогнутый; скорость фильтрации, как вертикальная, так и горизонтальная, выше.
Таким образом, можно сделать вывод, что при охлаждении ступенчато сужающихся пористых тепловыделяющих элементов возможны как установление стационарного режима, так и неограниченный разогрев, ведущий к расплавлению твердой фазы. Ступенчатое сужение элемента в верхней части приводит к существенному ухудшению охлаждения: основная часть газа проходит через среднюю часть элемента, где нет препятствий, поэтому ее сужение приводит к уменьшению суммарного расхода газа и, соответственно, к росту температуры и изменению других искомых величин. Наибольший разогрев твердой среды происходит в углах, образованных соединением нижних вертикальных стенок с горизонтальными непроницаемыми стенками.
Сравнивая охлаждение плавно и ступенчато сужающихся пористых тепловыделяющих элементов, можно заметить, что при ступенчатом сужении охлаждение хуже: в этом случае газ течет через элемент менее интенсивно, что влечет за собой рост температуры и изменение других искомых величин.
Заключение
Для моделирования двумерных нестационарных течений газа через твердые пористые тепловыделяющие элементы предложен численный метод, основанный на комбинации явных и неявных конечно-разностных схем. С применением этого метода решена задача о включении принудительной фильтрации в момент начала тепловыделения в твердой фазе для тепловыделяющих элементов плавно сужающейся и ступенчато сужающейся формы. Показано, что возможны как установление стационарного режима охлаждения, так и неограниченный разогрев пористого элемента, ведущий к расплавлению твердой фазы и нарушению описанного процесса. Также показано, что форма тепловыделяющего элемента сильно влияет на процесс его охлаждения: как плавное, так и ступенчатое сужение элемента в верхней части приводит к существенному ухудшению охлаждения. Обнаружены наиболее разогреваемые зоны у рассмотренных пористых элементов.
Список литературы
[1] Млслов В.П., Мясников В.П., Данилов В.Г. Математическое моделирование аварий-
ного блока Чернобыльской АЭС. М.: Наука, 1987. 144 с.
[2] Левин В.А., Луценко Н.А. Стационарный режим фильтрационного охлаждения пористого тепловыделяющего элемента // Юбилейный сб.: К тридцатилетию Института автоматики и процессов управления ДВО РАН. Владивосток: ИАПУ ДВО РАН, 2001. С. 151-159.
[3] Луценко Н.А. Одномерный стационарный режим фильтрации газа через слой неподвижного тепловыделяющего конденсированного материала // Дальневосточный мат. журн. 2002. Т. 3, № 1. С. 123-130.
[4] Луценко Н.А. Нестационарные режимы охлаждения пористого тепловыделяющего элемента // Мат. моделирование. 2005. Т. 17, № 3. С. 120-128.
[5] Левин В.А., Луценко Н.А. Возникновение неустойчивых режимов охлаждения пористого тепловыделяющего элемента при докритических краевых условиях // Горение и плазмохи-мия. 2005. Т. 3, № 2. С. 81-92.
[6] Левин В.А., Луценко Н.А. Течение газа через пористую тепловыделяющую среду при учете температурной зависимости вязкости газа // Инженерно-физ. журн. 2006. Т. 79, № 1. С. 35-40.
[7] НигмАтулин Р.И. Основы механики гетерогенных сред. М.: Наука, 1978. 336 с.
[8] Жумагулов Б.Т., Монахов В.Н. Гидродинамика нефтедобычи. Алматы: КазгосИНТИ, 2001. 336 с.
[9] Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидромеханика и теплообмен. В 2-х т. Т. 1: Пер. с англ. М.: Мир, 1990. 384 с.
[10] Годунов С.К., Рябенький В.С. Разностные схемы. Введение в теорию. М.: Наука, 1977. 439 с.
Поступила в редакцию 10 мая 2006 г., в переработанном виде — 28 июля 2006 г.
