Сохрани и опубликуйсвоё исследование
О проекте | Cоглашение | Партнёры
магистерская диссертация по направлению подготовки : 15.04.03 - Прикладная механика
Источник: Федеральное государственное автономное образовательное учреждение высшего образования «Дальневосточный федеральный университет»
Комментировать 0
Рецензировать 0
Скачать - 3,2 МБ
Enter the password to open this PDF file:
-
Оглавление Введение ............................................................................................................... 2 Глава 1. Основные уравнения математической модели ................................. 7 1.1. Постановка задачи .................................................................................... 7 1.2 Математическая модель ............................................................................ 8 1.3. Начальные условия ................................................................................. 11 1.4. Обезразмеривание уравнений ............................................................... 13 1.5. Численый метод ...................................................................................... 18 Глава 2. Результаты решения ........................................................................... 28 2.1. Тестовая задача ....................................................................................... 28 2.2. Результаты вычислительных экспериментов о нестационарном режиме нагревания пористого элемента ..................................................... 33 2.3. Результаты вычислительных экспериментов о нестационарном режиме охлаждения пористого элемента .................................................... 41 Заключение ........................................................................................................ 51 Литература ......................................................................................................... 52
Введение Сегодня быстрое развитие возобновляемых источников энергии (особенно солнечной и ветровой энергии) может эффективно смягчить серьезный энергетический кризис во всем мире. Однако из-за прерывистости и неравномерности производства возобновляемой энергии она не может полностью энергию. Накопитель заменить энергии традиционную сжатого воздуха тепловую (Воздушно- аккумулирующая газотурбинная электростанция, ВАГТЭ,Compressed Air Energy Storage, CAES) является одной из многообещающих систем хранения энергии, позволяющий максимально эффективно использовать энергетический потенциал электростанций, и может принести многочисленные выгоды для работы энергосистем, что в свою очередь, является хорошим дополнением к решению этой проблемы [13, 20]. Накопитель сжатого воздуха аккумулирует энергию в виде сжатого воздуха под высоким давлением в периоды низкой потребности электроэнергии, а затем высвобождает его, пропуская через турбину для выработки электроэнергии, чтобы удовлетворить высокий спрос в течение пикового периода времени. Это позволяет снизить стоимость электроэнергии за счет увеличения КПД, а также обеспечить её бесперебойную подачу потребителю. Первые из CAES были построены и введены в эксплуатацию в Германии (“Huntorf”, сдана в эксплуатацию в 1978 г.) и в США (Mclntosh, Алабама, 1991 г.). Обе установки в качестве хранилища используют подземные пещеры, полученные после растворения каменной соли через буровую скважину. Однако такие станции теряют большую часть потенциальной энергии в виде тепла, вырабатываемого на этапе заряда от сжатия воздуха. Затем на этапе разряда расширяющийся газ подогревают энергией от сжигания природного газа. Перспективным развитием в этой области, которое значительно уменьшает потери тепла, является добавление в систему теплообменника, аккумулирующего тепло 2
от сжатия, а затем использующего его при расширении воздуха [21], такое CAES называется адиабатическим (A-CAES). А за счет повторного использования накопленной тепловой энергии в A-CAES исключается загрязнение атмосферы от горения углеводородов. Рабочий процесс установки по хранению сжатого воздуха с добавлением теплообменника описывается следующим образом (рис. 1). Во время режима сжатия избыточное электричество, выработанное на электростанции, используется для управления цепью компрессоров. Воздух, проходя через них в процессе сжатия, выделяет тепловую энергию, после чего, проходя через теплообменник, отдает излишки тепловой энергии окружающей среде (стенкам пористой среды). После этого воздух попадает в хранилище сжатого воздуха, в качестве которого может быть подземная пещера или баллоны. Сжатый воздух хранится при высоком давлении и при температуре окружающего пласта. В режиме расширения сжатый воздух, сохраненный под высоким давлением, выпускается и нагревается при прохождении теплообменника (от тепла, накопленного при сжатии), и затем расширяется при прохождении группы турбин. Таким образом, энергия, хранящаяся в сжатом воздухе, преобразуется в электрическую без участия процессов горения. Теоретически общая эффективность адиабатического хранилища энергии сжатого воздуха выше, чем у традиционной технологии хранения сжатого воздуха, поскольку адиабатическим усовершенствованные сжатым воздухом системы повторно хранения используют с тепло, генерируемое в процессе сжатия [18,21]. 3
КОМПРЕССОР ТУРБИНА ХРАНИЛИЩЕ ТЕПЛОВОЙ ЭНЕРГИИ ХРАНИЛИЩЕ СЖАТОГО ВОЗДУХА РЕЖИМ СЖАТИЯ РЕЖИМ РАСШИРЕНИЯ Рис. 1 Схема накопителя энергии сжатого воздуха В связи с этим пристальный интерес вызывает анализ внутренних задач конвекции в теплообменнике, для выявления наиболее эффективного решения при его проектировании. Поскольку наполнитель теплообменника в CAES является пористой средой, то движение газа через него можно рассматривать как движение двухкомпонентной гетерогенной среды и описывать методами механики сплошных гетерогенных сред, тогда математическую модель процесса можно построить в рамках модели взаимодействующих континуумов [11,12]. Основы для исследования движения жидкостей через пористые среды, то есть тела, пронизанные системой сообщающихся между собой пустот, были заложены французским инженером Анри Дарси. Он в 1856 г. на основе опытов по фильтрации жидкости через цилиндр, заполненный песком, экспериментально установил, что скорость фильтрации жидкости при движении через пористую среду пропорциональна проницаемости среды и градиенту давления и обратно пропорционально вязкости жидкости. При этом скорость фильтрации направлена в сторону убывания 4
давления [2,9,19]. Согласно механике сплошных сред [16,17] между жидкостью и газом мало различий и одни и те же законы фильтрации справедливы для обоих состояний. Большое влияние в развитие теории фильтрации оказали работы Н.Е. Жуковского [3], Н.Н. Павловского [14]. Значительное влияние оказали работы П.Я. Полубариновой-Кочиной [15], чьи исследования в двух сборниках использовались для решения задач гидроэнергетики, нефте- и газодобычи. В работах [6,7], посвященных исследованию нестационарных процессов охлаждения газом пористых тепловыделяющих объектов, были предложены математическая модель и численный метод, основанный на комбинации явных и неявных конечно-разностных схем. Дальнейшие исследования получили развитие в направление моделирования движения газа через слой теплоаккумулирующего материала с фазовым переходом, где используется ранее описанный метод с адаптацией к новой модели, а полученные численные расчеты при сравнении с экспериментальными данными демонстрируют хорошее совпадение [4]. Цель настоящей работы состоит в том, чтобы получить и проанализировать решение задачи о нестационарном охлаждении и нагреве газа пористым объектом. Для достижения цели поставлены и решены следующие задачи: 1) Изучение литературы, связанной с темой работы. 2) Выбор математической модели для исследования нестационарного течения газа через пористые среды в вертикальном цилиндре. 3) Выбор или разработка численного метода для решения уравнений модели. 4) Разработка программного средства, реализующего предложенный численный метод, и его тестирование. 5) Получение численного решения поставленной задачи и проведение анализа полученных результатов. 5
Методы исследования Математическая модель базируется на классических подходах механики сплошных многокомпонентных сред основанных на идее взаимодействующих континуумов. Численный метод основан на методе конечных разностей с комбинацией явных и неявных конечно-разностных схем, с применением демпфирующих членов для уменьшения дисперсионной ошибки. Для реализации численного метода были написаны алгоритмы на языке С++ в программном продукте Microsoft Visual Studio и MatLab. 6
Глава 1. Основные уравнения математической модели 1.1. Постановка задачи Объектом исследования является неподвижный, однородный, пористый объект высотой Н, который имеет боковые непроницаемые стенки и пористые элемент, сквозь который может проходить газ. Вертикальные и ровные боковые стенки закрывают пористый материал от окружающей среды и исключают теплопроводность в пристеночных областях. Снизу и сверху объект открыт и в его нижнюю часть с заданным давлением и температурой подается горячий газ (или холодный), который в процессе движения через пористую среду нагревает (или охлаждает) её в результате теплообмена и выходит охлажденным (или нагретым) в открытое пространство (рис. 1.1.1). Рис. 1.1.1 Пористый теплоаккумулирующий объект: a) процесс аккумулирования тепла, б) процесс разряда теплоаккумулятора 7
Пологая, что по горизонтальной переменной отсутствуют изменения параметров и пренебрегая эффектами в пристеночных областях, можно рассматривать задачу относительно одной пространственной переменной изменяющейся по вертикале. 1.2 Математическая модель Математическую модель описанного процесса будем рассматривать в качестве движения двухкомпонентной сплошной среды, состоящей из твердой и газовой фазы, в рамках модели взаимопроникающего движения двух взаимодействующих континуумов, которые представляют каждую из фаз. Это означает, что при рассмотрении некоторого объема пористой среды, который содержит в себе частицы твердой среды и газа, мы мысленно “размазываем” эти две среды по всему рассматриваемому объёму и полагаем, что в каждой точке этого объема содержится как твердая фаза, так и газовая. Данный подход излагается в работах Р.И. Нигматулина [11-12] и используется в механике сплошных гетерогенных сред, а классическим примером служит теория фильтрации. Так как фазовые переходы отсутствуют, то масса каждой из фаз сохраняется, и закон сохранения массы для газа можно записать в следующем виде: 𝑎 𝜕𝜌𝑔 𝜕𝑎𝜌𝑔 𝑣𝑔 + =0 𝜕𝑡 𝜕𝑥 (1.2.1) где 𝑡 − время, с; 𝑥 - пространственная переменная, м; 𝜌𝑔 – истинная плотность газовой фазы, кг/м3 ; 𝜐𝑔 - скорость движения газовой фазы, м/с; 𝑎 - пористость (относительный объём газовой фазы в двухфазной смеси). 8
Относительный объём, занятый твердой фазой, будет равен (1 – a). С учетом течения газ в поле силы тяжести уравнение движения принимает следующий вид: 𝜕𝑣𝑔 𝜕𝑣𝑔 𝜕𝑝 𝜇 𝜌𝑔 (1 + 𝜒(1 − 𝑎)) ( + 𝑣𝑔 − 𝜌𝑔 𝑔 − 𝑎 𝑣𝑔 )=− 𝜕𝑡 𝜕𝑥 𝜕𝑥 𝑘1 (1.2.2) где 𝑝 – давление газа, Па; 𝜇 – коэффициент динамической вязкости газа, Па ∙ с; 𝑔 – ускорение силы тяжести, м2 /2; 𝑘1 – коэффициент проницаемости твердой фазы, м2 ; 𝜒 – коэффициент присоединенной массы, учитывающий инерционное взаимодействие фаз при их ускоренном относительном движении. Уравнение энергии твердой фазы с учётом того, что она неподвижна, а конвективный теплообмен между фазами пропорционален разности температур, имеет следующий вид: 𝜕𝑇𝑐 𝜕 2 𝑇𝑐 (1 − 𝑎)𝜌𝑐 𝐶𝑐 = −𝛼(𝑇𝑐 − 𝑇𝑔 ) + (1 − 𝑎)𝜆𝑐 𝜕𝑡 𝜕𝑥 2 (1.2.3) где 𝑇𝑔 – температура газа, К; 𝑇𝑐 – температура твердой фазы, К; 𝐶𝑐 – теплоёмкость твердой фазы, Дж/(кг ∙ К); 𝜆𝑐 – коэффициент теплопроводности твердой фазы, Вт/(м ∙ К); 𝜌𝑐 – плотность конденсированной фазы, кг/м3 ; 𝛼 – константа, определяющая интенсивность межфазового теплообмена, Вт(м3 ∙ К). Уравнение энергии газа: 𝜕𝑇𝑔 𝜕𝑇𝑔 𝜕𝑝 𝜕𝑝 𝜇 𝑎𝜌𝑔 𝐶𝑔𝑝 ( + 𝑣𝑔 ) = 𝛼(𝑇𝑐 − 𝑇𝑔 ) + 𝑎( + 𝑣𝑔 ) + 𝑎2 𝑣𝑔 2 𝜕𝑡 𝜕𝑥 𝜕𝑡 𝜕𝑥 𝑘1 (1.2.4) 9
где 𝐶𝑔𝑝 – теплоёмкость газа, Дж/(кг ∙ К). Для замыкания системы уравнений, описывающих движение газа в пористой среде, добавим уравнение состояния совершенного газа: 𝑝 = 𝜌𝑔 𝑅𝑇𝑔 (1.2.5) где 𝑅 – газовая постоянная, Дж/(кг ∙ К). Для твердой фазы уравнения движения и сохранения массы вырождаются, поскольку она неподвижна и имеет постоянную плотность. Таким образом, конечная система уравнений состоит из пяти уравнений (1.2.1 − 1.2.5) с пятью неизвестными: 𝜕𝑇𝑐 𝜕 2 𝑇𝑐 (1 − 𝑎)𝜌𝑐 𝐶𝑐 = −𝛼(𝑇𝑐 − 𝑇𝑔 ) + (1 − 𝑎)𝜆𝑐 𝜕𝑡 𝜕𝑥 2 𝜕𝑇𝑔 𝜕𝑇𝑔 𝜕𝑝 𝜕𝑝 𝜇 𝑎𝜌𝑔 𝐶𝑔𝑝 ( + 𝑣𝑔 ) = 𝛼(𝑇𝑐 − 𝑇𝑔 ) + 𝑎( + 𝑣𝑔 ) + 𝑎2 𝑣𝑔 2 𝜕𝑡 𝜕𝑥 𝜕𝑡 𝜕𝑥 𝑘1 𝜕𝑣𝑔 𝜕𝑣𝑔 𝜕𝑝 𝜇 𝜌𝑔 (1 + 𝜒(1 − 𝑎)) ( + 𝑣𝑔 − 𝜌𝑔 𝑔 − 𝑎 𝑣𝑔 )=− 𝜕𝑡 𝜕𝑥 𝜕𝑥 𝑘1 𝜕𝜌𝑔 𝜕𝑎𝜌𝑔 𝑣𝑔 𝑎 + =0 𝜕𝑡 𝜕𝑥 R 𝑝 = 𝜌 𝑇 𝑔 { M 𝑔 При исследовании динамическая вязкость, неизотермической как правило, фильтрации полагается (1.2.6) жидкости зависящей от температуры, однако при моделировании неизотермической фильтрации газа вязкость часто принимается постоянной. Далее будем считать, что динамическая вязкость газа зависит от температуры по формуле Сазерленда [5]: 𝑇𝑔1.5 𝜇 = с𝑠1 . 𝑐𝑠2 + 𝑇𝑔 (1.2.7) 10
На открытых границах (входе и выходе) пористого объекта известны условия теплообмена и значение давления. Также на входе известна температура поступающего газа. Таким образом, краевые условия для системы уравнений (2.1.6) можно записать следующим образом: 𝑝(𝑥)|𝑥=0 = 𝑝0 , 𝑝(𝑥)|𝑥=𝐻 = 𝑝𝐻 , 𝑇𝑔 (𝑥)| = 𝑇𝑔0 , 𝑥=0 𝑑𝑇𝑐 (𝑥) 𝛼∗ = (𝑇𝑐0 − 𝑇𝑔0 ), | 𝑑𝑥 𝑥=0 𝜆𝑐 𝑑𝑇𝑐 (𝑥) = 0, | 𝑑𝑥 𝑥=𝐻 (1.2.8) где 𝑝0 , 𝑝𝐻 – известные величины давления на границах пористого элемента; 𝑇𝑔0 – температура газа поступающего в нижнюю границу пористого элемента; 𝛼∗ - коэффициент теплоотдачи, Дж/(м2 ∙ К ∙ с). 1.3. Начальные условия Для получения расчетных данных для начального момента времени, рассмотрим стационарную задачу с условиями, что в начальный момент никаких процессов не происходит, давление на входе в объект и на выходе из него соответствует атмосферному давлению на заданных высотах, следовательно, движение воздуха в объекте отсутствует, поэтому скорость газа, а следовательно и скорость фильтрации газа 𝑣𝑔 = 0. Температура газа и теплообменника равна температуре окружающей среды. Тогда система (1.2.6) принимает следующий вид: 0 = −𝛼(𝑇𝑐 − 𝑇𝑔 ) 𝜕𝑝 = −𝜌𝑔 𝑔 𝑢=0 𝜕𝑥 𝜕𝑝 ⟹ 𝑀𝑝 . 0=− − 𝜌𝑔 𝑔 𝜌𝑔 = 𝜕𝑥 𝑅𝑇𝑔 { { 𝑝 = 𝜌𝑔 𝑅𝑇𝑔 11
Подставляем формулу для плотности в дифференциальное уравнение для 𝜕𝑝 𝜕𝑥 : 𝜕𝑝 𝑀𝑝𝑔 =− . 𝜕𝑥 𝑅𝑇𝑔 В результате мы получаем дифференциальное уравнение, описывающее давление газа 𝑝 как функцию переменной 𝑥. 𝜕𝑝 𝑀𝑔 =− 𝜕𝑥. 𝑝 𝑅𝑇𝑔 Интегрирование приводит к следующему уравнению: ln |𝑝| = − 𝑀𝑔 𝑥 + ln С. 𝑅𝑇𝑔 Избавляясь от логарифмов, получаем уравнение, из которого следует, что в начальный момент во всей системе давление распределено согласно барометрической формуле. 𝑝 = С𝑒 Константа C определяется 𝑀𝑔𝑥 − 𝑅𝑇𝑔 из . начального условия 𝑝(𝑥)|𝑥=0 = 𝑝0 , где 𝑝0 – это среднее атмосферное давление над уровнем моря. Таким образом, зависимость атмосферного давления газа от высоты выражается формулой: 𝑝 = 𝑝0 𝑒 𝑀𝑔𝑥 − 𝑅𝑇𝑔 , а система уравнений, описывающая состояние в начальный момент, принимает вид: 𝑇𝑐 = 𝑇𝑔 𝑢=0 𝑝 = 𝑝0 𝑒 { 𝜌𝑔 = 𝑀𝑔𝑥 − 𝑅𝑇𝑔 . (1.3.1) 𝑝 𝑅𝑇𝑔 12
1.4. Обезразмеривание уравнений Рассмотрим нестационарное одномерное движение газа. Будем решать систему (1.2.7) c граничными условиями (1.2.8). Введём понятие скорости фильтрации газа 𝑢 = 𝑎𝜐𝑔 , которое в дальнейшем будем использовать вместо скорости газа, и перепишем систему уравнений движения газа через пористую среду с её учетом: 𝜕𝑇𝑐 𝜕 2 𝑇𝑐 = −𝛼(𝑇𝑐 − 𝑇𝑔 ) + (1 − 𝑎)𝜆𝑐 𝜕𝑡 𝜕𝑥 2 𝜕𝑇𝑔 𝜕𝑇𝑔 𝜕𝑝 𝜕𝑝 𝜇 2 𝜌𝑔 𝐶𝑔𝑝 (𝑎 +𝑢 +𝑢 + 𝑢 ) = 𝛼(𝑇𝑐 − 𝑇𝑔 ) + 𝑎 𝜕𝑡 𝜕𝑥 𝜕𝑡 𝜕𝑥 𝑘1 (1 − 𝑎)𝜌𝑐 𝐶𝑐 { 𝜌𝑔 (1 + 𝜒(1 − 𝑎)) 𝜕𝑢 𝑢 𝜕𝑢 𝜕𝑝 𝜇 − 𝜌𝑔 𝑔 − 𝑢 ( + )=− 𝑎 𝜕𝑡 𝑎 𝜕𝑥 𝜕𝑥 𝑘1 𝜕𝜌𝑔 𝜕𝜌𝑔 𝑢 𝑎 + =0 𝜕𝑡 𝜕𝑥 R 𝑝 = 𝜌𝑔 𝑇𝑔 M (1.4.1) Прежде чем перейти к созданию алгоритма для решения уравнений математической модели (1.3.1), необходимо привести эти уравнения к безразмерному виду, т.е. провести операцию обезразмеривания переменных, в результате которой все переменные математической модели будут иметь одинаковый порядок. Для этого введем безразмерные переменные через соотношения: 𝑡 = 𝑡∗ 𝑡̃, 𝑥 = 𝐻𝑥̃, 𝑢 = 𝑢∗ 𝑢̃, 𝑝 = 𝑝∗ 𝑝̃, (1.4.2) 𝜌𝑔 = 𝜌∗ 𝜌̃𝑔 , 𝑇𝑐 = 𝑇∗ 𝑇̃𝑐 , 13
𝑇𝑔 = 𝑇∗ 𝑇̃𝑔 , где 𝑡∗ , 𝑢∗ − характерные значения времени и скорости фильтрации газа; 𝑝∗ , 𝜌∗ , 𝑇∗ − давление, плотность, температура газа при «нормальных» условиях, а знак тильда означает безразмерную величину. Подставим (1.4.2) в исходную систему уравнений (1.4.1) и условия (1.2.8) и получим систему уравнений: 𝜕𝑇∗ 𝑇̃𝑐 𝜕 2 𝑇∗ 𝑇̃𝑐 ̃ ̃ (1 = −𝛼(𝑇∗ 𝑇𝑐 − 𝑇∗ 𝑇𝑔 ) + − 𝑎)𝜆𝑐 𝜕𝑡∗ 𝑡̃ 𝜕(𝐻𝑥̃)2 𝜕𝑇∗ 𝑇̃𝑔 𝜕𝑇∗ 𝑇̃𝑔 𝑢∗2 с𝑠1 √𝑇∗ 𝑇̃𝑔1.5 ̃ ̃ 𝜌∗ 𝜌̃𝑔 𝐶𝑔𝑝 (𝑎 + 𝑢∗ 𝑢̃ 𝑢̃2 ) = 𝛼(𝑇∗ 𝑇𝑐 − 𝑇∗ 𝑇𝑔 ) + 𝜕𝑡∗ 𝑡̃ 𝜕𝐻𝑥̃ 𝑘1 𝑐̃𝑠2 + 𝑇̃𝑔 (1 − 𝑎)𝜌𝑐 𝐶𝑐 𝜌∗ 𝜌̃𝑔 (1 + 𝜒(1 − 𝑎)) 𝜕𝑢∗ 𝑢̃ 𝑢∗ 𝑢̃ 𝜕𝑢∗ 𝑢̃ 𝜕𝑝∗ 𝑝̃ 𝜇 + − 𝜌∗ 𝜌̃𝑔 𝑔 − 𝑢∗ 𝑢̃ ( )=− 𝑎 𝜕𝑡∗ 𝑡̃ 𝑎 𝜕𝐻𝑥̃ 𝜕𝐻𝑥̃ 𝑘1 𝜕𝜌∗ 𝜌̃𝑔 𝜕𝜌∗ 𝜌̃𝑔 𝑢∗ 𝑢̃ 𝑎 + =0 𝜕𝑡∗ 𝑡̃ 𝜕𝐻𝑥̃ 𝑝∗ 𝑝̃ = 𝜌∗ 𝜌̃𝑔 R𝑇∗ 𝑇̃𝑔 { (1.4.3) и граничные условия: 𝑝̃(𝑥)|𝑥=0 = 𝑝̃0 , 𝑝̃(𝑥)|𝑥=𝐻 = 𝑝̃𝐻 , 𝑇̃𝑔 (𝑥)| = 𝑇̃𝑔0 , 𝑥=0 𝑑𝑇̃с (𝑥) 𝛼∗ 𝐻 = | (𝑇̃с |𝑥=0 − 𝑇̃𝑔0 ), 𝑑𝑥̃ 𝑥=0 𝜆𝑐 𝑑𝑇̃с (𝑥) (𝑥)| = 0, 𝑑𝑥̃ 𝑥=𝐻 (1.4.4) После несложных преобразований приведем (1.4.3) к виду: (1 − 𝑎)𝜌𝑐 𝐶𝑐 𝑇∗ 𝜕𝑇̃𝑐 𝑇∗ 𝜕 2 𝑇̃𝑐 ̃ ̃ (1 = −𝛼𝑇∗ (𝑇𝑐 − 𝑇𝑔 ) + − 𝑎)𝜆𝑐 2 𝑡∗ 𝜕𝑡̃ 𝐻 𝜕𝑥̃ 2 14
𝑇∗ 𝜕𝑇̃𝑔 𝑇∗ 𝜕𝑇̃𝑔 𝑎 𝜕𝑝̃ 𝜌∗ 𝜌̃𝑔 𝐶𝑔𝑝 (𝑎 + 𝑢∗ 𝑢̃ + ) = 𝛼𝑇∗ (𝑇̃𝑐 − 𝑇̃𝑔 ) + 𝑡∗ 𝜕𝑡̃ 𝐻 𝜕𝑥̃ 𝑝∗ 𝑡∗ 𝜕𝑡̃ 𝑢∗2 с𝑠1 √𝑇∗ 𝑇̃𝑔1.5 𝑢∗ 𝑝∗ 𝜕𝑝̃ + 𝑢̃ ++ 𝑢̃2 𝐻 𝜕𝑥̃ 𝑘1 𝑐̃𝑠2 + 𝑇̃𝑔 (1.4.5) 𝜌∗ 𝜌̃𝑔 (1 + 𝜒(1 − 𝑎)) 𝑢∗ 𝜕𝑢̃ 𝑢∗ 𝑢̃ 𝑢∗ 𝜕𝑢̃ 𝜕𝑝∗ 𝑝̃ 𝜇 + − 𝜌∗ 𝜌̃𝑔 𝑔 − 𝑢∗ 𝑢̃ ( )=− 𝑎 𝑡∗ 𝜕𝑡̃ 𝑎 𝐻 𝜕𝑥̃ 𝜕𝐻𝑥̃ 𝑘1 𝑎 𝜌∗ 𝜕𝜌̃𝑔 𝜌∗ 𝑢∗ 𝜕𝜌̃𝑔 𝑢̃ + =0 𝑡∗ 𝜕𝑡̃ 𝐻 𝜕𝑥̃ 𝑝̃ = 𝜌̃𝑔 𝑇̃𝑔 Рассмотрим первое уравнение (1.4.5) и поделим его на (1 − 𝑎)𝜌𝑐 𝐶𝑐 : 𝜕𝑇̃𝑐 𝛼𝑡∗ 𝜆𝑐 𝑡∗ 𝜕 2 𝑇̃𝑐 =− . (𝑇̃ − 𝑇̃𝑔 ) + (1 − 𝑎)𝜌𝑐 𝐶𝑐 𝑐 𝜕𝑡̃ 𝜌𝑐 𝐶𝑐 𝐻2 𝜕𝑥̃ 2 (1.4.6) Запишем произведения констант в качестве параметров подобия: 𝜋1 = 𝛼𝑡∗ , (1 − 𝑎)𝜌𝑐 𝐶𝑐 (1.4.7) 𝜆𝑐 𝑡∗ . 𝜌𝑐 𝐶𝑐 𝐻2 (1.4.8) 𝜋2 = Рассмотрим второе уравнение (1.4.5) и раскроем скобки в левой части: 𝑇∗ 𝜕𝑇̃𝑔 𝑇∗ 𝜕𝑇̃𝑔 𝜌∗ 𝜌̃𝑔 𝐶𝑔𝑝 𝑎 + 𝜌∗ 𝜌̃𝑔 𝐶𝑔𝑝 𝑢∗ 𝑢̃ 𝑡∗ 𝜕𝑡̃ 𝐻 𝜕𝑥̃ 𝑎 𝜕𝑝̃ 𝑢∗ 𝑝∗ 𝜕𝑝̃ 𝑢∗2 с𝑠1 √𝑇∗ 𝑇̃𝑔1.5 = 𝛼𝑇∗ (𝑇̃𝑐 − 𝑇̃𝑔 ) + + 𝑢̃ + 𝑢2 . ̃ ̃ 𝑝∗ 𝑡∗ 𝜕𝑡 𝐻 𝜕𝑥̃ 𝑘1 𝑐̃𝑠2 + 𝑇𝑔 Перенесем 𝜌∗ 𝜌̃𝑔 𝐶𝑔𝑝 𝑢∗ 𝑢̃ 𝑇∗ 𝜕𝑇̃𝑔 𝐻 𝜕𝑥̃ в правую часть и поделим всё уравнение на 𝜌∗ 𝜌̃𝑔 𝐶𝑔𝑝 𝑎: 15
𝜕𝑇̃𝑔 𝑢∗ 𝑡∗ 𝜕𝑇̃𝑔 𝛼𝑡∗ (𝑇̃𝑐 − 𝑇̃𝑔 ) 𝑝∗ 𝜕𝑝̃ =− 𝑢̃ + + 𝜕𝑡̃ 𝑎𝐻 𝜕𝑥̃ 𝜌∗ 𝐶𝑔𝑝 𝑎 𝜌̃𝑔 𝜌∗ 𝐶𝑔𝑝 𝑇∗ 𝜌̃𝑔 𝜕𝑡̃ 𝑇̃𝑔1.5 𝑢2 𝑡∗ 𝑢∗ 𝑝∗ 𝜕𝑝̃ 𝑢̃ с𝑠1 𝑢∗2 𝑡∗ + + 𝜌∗ 𝐶𝑔𝑝 𝑎𝑇∗ 𝐻 𝜕𝑥̃ 𝜌̃𝑔 𝑘1 𝑎𝜌∗ 𝐶𝑔𝑝 √𝑇∗ 𝑐̃𝑠2 + 𝑇̃𝑔 𝜌̃𝑔 (1.4.9) Запишем произведения констант в качестве параметров подобия: 𝜕𝑇̃𝑔 𝜕𝑇̃𝑔 𝑇̃𝑔1.5 𝑢2 (𝑇̃𝑐 − 𝑇̃𝑔 ) 𝜕𝑝̃ 𝜕𝑝̃ 𝑢̃ = −𝜋3 𝑢̃ + 𝜋4 + 𝜋5 + 𝜋3 𝜋5 + 𝜋6 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝜌̃𝑔 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝑐̃𝑠2 + 𝑇̃𝑔 𝜌̃𝑔 𝜋3 = 𝑢∗ 𝑡∗ , 𝑎𝐻 (1.4.10) 𝜋4 = 𝛼𝑡∗ , 𝑎𝜌∗ 𝐶𝑔𝑝 (1.4.11) 𝑝∗ , 𝜌∗ 𝐶𝑔𝑝 𝑇∗ (1.4.12) 𝜋5 = 𝜋6 = с𝑠1 𝑢2∗ 𝑡∗ 𝑘1 𝑎𝜌∗ 𝐶𝑔𝑝 √𝑇∗ . (1.4.13) Рассмотрим третье уравнение (1.4.5) и раскроем скобки в левой части: 𝜌∗ 𝜌̃𝑔 (1 + 𝜒(1 − 𝑎)) 𝑢∗ 𝜕𝑢̃ 𝜌∗ 𝜌̃𝑔 (1 + 𝜒(1 − 𝑎)) 𝑢∗ 𝑢̃ 𝑢∗ 𝜕𝑢̃ + 𝑎 𝑡∗ 𝜕𝑡̃ 𝑎 𝑎 𝐻 𝜕𝑥̃ с𝑠1 𝑢∗ √𝑇∗ 𝑇̃𝑔1.5 𝜕𝑝∗ 𝑝̃ =− − 𝜌∗ 𝜌̃𝑔 𝑔 − 𝑢̃ 𝜕𝐻𝑥̃ 𝑘1 𝑐𝑠2 + 𝑇̃𝑔 Перенесем уравнение на ̃𝑔 (1+𝜒(1−𝑎)) 𝑢∗ 𝑢 𝜌∗ 𝜌 ̃ 𝑢∗ 𝜕𝑢 ̃ 𝑎 𝐻 𝜕𝑥̃ 𝑎 в правую часть и поделим всё ̃𝑔 (1+𝜒(1−𝑎)) 𝑢∗ 𝜌∗ 𝜌 𝑎 𝑡∗ 16
𝜕𝑢̃ 𝜕𝑢̃ 𝑎𝑡∗ 𝑝∗ 𝜕𝑝̃ 𝑎𝑡∗ 𝑔 = −𝜋3 𝑢̃ − − 𝜕𝑡̃ 𝜕𝑥̃ 𝜌∗ 𝑢∗ 𝐻(1 + 𝜒(1 − 𝑎))𝜌̃𝑔 𝜕𝑥̃ 𝑢∗ (1 + 𝜒(1 − 𝑎)) 𝑢̃ 𝑇̃𝑔1.5 − . 𝑘1 𝜌∗ (1 + 𝜒(1 − 𝑎)) 𝜌̃𝑔 𝑐𝑠2 + 𝑇̃𝑔 𝑎𝑡∗ с𝑠1 √𝑇∗ (1.4.14) Запишем произведения констант в качестве параметров подобия: 𝜕𝑢̃ 𝜕𝑢̃ 𝜋7 𝜕𝑝̃ 𝑢̃ 𝑇̃𝑔1.5 = −𝜋3 𝑢̃ − − 𝜋8 − 𝜋9 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝜕𝑥̃ 𝜌̃𝑔 𝑐𝑠2 + 𝑇̃𝑔 𝜋7 = 𝑝∗ 𝑎𝑡∗ 𝐻𝜌∗ 𝑢∗ (1 + 𝜒(1 − 𝑎)) 𝜋8 = 𝜋9 = 𝑔𝑎𝑡∗ 𝑢∗ (1 + 𝜒(1 − 𝑎)) , , 𝑎𝑡∗ с𝑠1 √𝑇∗ 𝑘1 𝜌∗ (1 + 𝜒(1 − 𝑎)) (1.4.15) (1.4.16) . (1.4.17) Рассмотрим третье уравнение системы (1.4.5) и переведем плотность газа согласно уравнению Клапейрона-Менделеева: 𝑎𝑝∗ 𝜕 𝑝̃ 𝑝∗ 𝑢∗ 𝜕 𝑝̃𝑢̃ ( )+ ( ) = 0. 𝑅𝑇∗ 𝑡∗ 𝜕𝑡̃ 𝑇̃𝑔 𝑅𝑇∗ 𝐻 𝜕𝑥̃ 𝑇̃𝑔 Перенесем 𝑎𝑝∗ 𝑅𝑇∗ 𝑡∗ ̃ 𝑝∗ 𝑢∗ 𝜕 𝑝̃𝑢 ( ) в правую часть и поделим всё уравнение на 𝑅𝑇∗ 𝐻 𝜕𝑥̃ 𝑇̃𝑔 : 𝜕 𝑝̃ 𝑢∗ 𝑡∗ 𝜕 𝑝̃𝑢̃ ( )=− ( ). 𝜕𝑡̃ 𝑇̃𝑔 𝑎𝐻 𝜕𝑥̃ 𝑇̃𝑔 (1.4.18) Таким образом, систему уравнений (1.4.5) с учетом преобразований (1.4.6 – 1.4.18) можно переписать в безразмерных переменных: 17
𝜕𝑇̃𝑐 𝜕 2 𝑇̃𝑐 ̃ ̃ = −𝜋1 (𝑇𝑐 − 𝑇𝑔 ) + 𝜋2 𝜕𝑡̃ 𝜕𝑥̃ 2 𝜕𝑇̃𝑔 𝜕𝑇̃𝑔 𝜋4 𝑇̃𝑔1.5 𝑢2 𝜕𝑝̃ 𝜕𝑝̃ 𝑢̃ ̃ ̃ = −𝜋3 𝑢̃ + (𝑇 − 𝑇𝑔 ) + 𝜋5 ( + 𝜋3 ) + 𝜋5 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝑐 𝜌̃𝑔 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝑐̃𝑠2 + 𝑇̃𝑔 𝜌̃𝑔 (1.4.19) 𝜕𝑢̃ 𝜕𝑢̃ 𝜋7 𝜕𝑝̃ 𝑢̃ = −𝜋3 𝑢̃ − − 𝜋8 − 𝜋9 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝜕𝑥̃ 𝜌̃𝑔 𝜕 𝜕 𝜌̃𝑔 = −𝜋3 𝜌̃ 𝑢̃ 𝜕𝑡̃ 𝜕𝑥̃ 𝑔 𝑝̃ = 𝜌̃𝑔 𝑇̃𝑔 { 1.5. Численый метод Система уравнений (1.4.19) с краевыми условиями (1.4.4), моделирующая нестационарные одномерные течения газа через пористый объект, является гиперболически-параболической системой уравнений и в общем случае не может быть решена аналитически, поэтому решение получим методом конечных разностей на равномерной сетке с шагом по пространству h и шагом по времени τ, равным 𝑟ℎ2 , где 𝑟 = 𝑐𝑜𝑛𝑠𝑡. Нижний индекс при искомых сеточных функциях будет обозначать продвижение по пространству, верхний индекс – продвижение по времени. Тогда можно записать систему конечноразностных уравнений, аппроксимирующую исходную систему (1.4.19) со вторым порядком точности по h и первым порядком по τ, и устойчивую при некотором ограничении на r. Уравнение энергии для твердой фазы в безразмерных переменных 𝜕𝑇̃𝑐 𝜕 2 𝑇̃𝑐 ̃ ̃ = −𝜋1 (𝑇𝑐 − 𝑇𝑔 ) + 𝜋2 ; 𝜕𝑡̃ 𝜕𝑥̃ 2 Преобразуем уравнение энергии твердой фазы в явное конечноразностное уравнение, из которого определим температуру твердой среды. Заменяем первую производную по времени и вторую производную по пространству конечными разностями: 𝑛+1 𝑛 𝜕𝑇̃𝑐 𝑇̃𝑐𝑚 − 𝑇̃𝑐𝑚 = ; 𝜕𝑡̃ 𝜏 18
𝑛 𝑛 𝑛 𝜕 2 𝑇̃𝑐 𝑇̃𝑐𝑚+1 − 2𝑇̃𝑐𝑚 + 𝑇̃𝑐𝑚−1 = ; 𝜕𝑥̃ 2 ℎ2 при этом получим систему линейных алгебраических уравнений: 𝑛 𝑛 𝑛+1 𝑛 𝑛 𝑇̃𝑐𝑚 − 𝑇̃𝑐𝑚 𝑇̃𝑐𝑚+1 − 2𝑇̃𝑐𝑚 + 𝑇̃𝑐𝑚−1 𝑛 𝑛 = −𝜋1 (𝑇̃𝑐𝑚 − 𝑇̃𝑔𝑚 ) + 𝜋2 . 𝜏 ℎ2 𝑛+1 Выведем 𝑇̃𝑐𝑚 и для избавления от знаменателя записи распишем 𝜏 как 𝑟ℎ2 : 𝑛 𝑛 𝑛+1 𝑛 (1 𝑛 𝑛 𝑇̃𝑐𝑚 = 𝑇̃𝑐𝑚 − 𝜋1 𝜏) + 𝜋1 𝜏𝑇̃𝑔𝑚 + 𝜋2 𝑟(𝑇̃𝑐𝑚+1 − 2𝑇̃𝑐𝑚 + 𝑇̃𝑐𝑚−1 ). (1.5.1) Уравнение энергии для газа в безразмерных переменных 𝜕𝑇̃𝑔 𝜕𝑇̃𝑔 𝜋4 𝑇̃𝑔1.5 𝑢2 ̃ ̃ = −𝜋3 𝑢̃ + (𝑇 − 𝑇𝑔 ) + 𝜋6 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝑐 𝑐̃𝑠2 + 𝑇̃𝑔 𝜌̃𝑔 Преобразуем уравнение энергии газовой фазы в явное конечноразностное уравнение, из которого определим температуру газа. Заменяем первую производную по времени конечными разностями, для температуры газа будем использовать трехточечную конечно-разностную аппроксимацию первой производной: 𝑛+1 𝑛 𝜕𝑇̃𝑔 𝑇̃𝑔𝑚 − 𝑇̃𝑔𝑚 = ; 𝜕𝑡̃ 𝜏 𝑛 𝑛−1 𝜕𝑝̃ 𝑝̃𝑚 − 𝑝̃𝑚 = ; 𝜕𝑡̃ 𝜏 𝑛 𝑛 𝜕𝑝̃ 𝑝̃𝑚+1 − 𝑝̃𝑚−1 = ; 𝜕𝑥̃ 2ℎ 𝑛 𝑛 𝑛 ∓ 4𝑇𝑔𝑚∓1 ± 𝑇𝑔𝑚∓2 𝜕𝑇̃𝑔 ±3𝑇𝑔𝑚 = ; 𝜕𝑥̃ 2ℎ при этом получим систему алгебраических уравнений: 19
𝑛 𝑛 𝑛 𝑛+1 𝑛 ±3𝑇𝑔𝑚 ∓ 4𝑇𝑔𝑚∓1 ± 𝑇𝑔𝑚∓2 𝑇̃𝑔𝑚 − 𝑇̃𝑔𝑚 𝜋4 𝑛 𝑛 = −𝜋3 𝑢̃ + 𝑛 (𝑇̃𝑐𝑚 − 𝑇̃𝑔𝑚 ) 𝜏 2ℎ 𝜌̃𝑔𝑚 𝑛 1.5 𝑛 𝑛 𝑛 𝑛−1 𝑛 𝑛 )2 (𝑇̃𝑔𝑚 ) (𝑢̃𝑚 𝑝̃𝑚 − 𝑝̃𝑚 𝑝̃𝑚+1 − 𝑝̃𝑚−1 𝑢̃𝑚 + 𝜋5 ( + 𝜋3 . 𝑛 𝑛 ) + 𝜋5 𝑛 𝑛 𝜌̃𝑔𝑚 𝜏 2ℎ 𝜌̃𝑔𝑚 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 𝑛+1 Оставим слева 𝑇̃𝑔𝑚 , а всё остальное перенесем вправо: 𝑛+1 𝑇̃𝑔𝑚 = 𝑛 𝑇̃𝑔𝑚 − 𝑛 𝜏𝜋3 𝑢̃𝑚 𝑛 𝑛 𝑛 ±3𝑇𝑔𝑚 ∓ 4𝑇𝑔𝑚∓1 ± 𝑇𝑔𝑚∓2 2ℎ + 𝜏𝜋4 𝑛 ̃ ̃𝑛 𝑛 (𝑇𝑐𝑚 − 𝑇𝑔𝑚 ) + 𝜌̃𝑔𝑚 𝑛 1.5 𝑛 𝑛−1 𝑛 𝑛 )2 ̃𝑔𝑚 (𝑇 ) (𝑢̃𝑚 𝑝̃𝑚 − 𝑝̃𝑚 𝜋5 𝜋3 𝜏 𝑛 𝑢 ̃ 𝑚 𝑛 +𝜋5 + (𝑝̃𝑚+1 − 𝑝̃𝑚−1 ) 𝑛 + 𝜏𝜋6 ; 𝑛 𝑛 𝑛 𝜌̃𝑔𝑚 2ℎ 𝜌̃𝑔𝑚 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 𝑛+1 𝑛 𝑇̃𝑔𝑚 = 𝑇̃𝑔𝑚 (1 − 𝜋4 𝜏 𝜋3 𝑟ℎ 𝑛 𝜋4 𝜏 𝑛 𝑛 𝑛 𝑛 ̃ − 𝑢 ̃ ∓ 4𝑇 ± 𝑇 + ) (±3𝑇 ) 𝑚 𝑔𝑚 𝑔𝑚∓1 𝑔𝑚∓2 𝑛 𝑛 𝑇𝑐𝑚 + 𝜌̃𝑔𝑚 2 𝜌̃𝑔𝑚 𝑛 1.5 𝑛 𝑛−1 𝑛 𝑛 )2 ̃ (𝑇 ) (𝑢̃𝑚 𝑝̃𝑚 − 𝑝̃𝑚 𝜋5 𝜋3 𝑟ℎ 𝑛 𝑢 ̃ 𝑔𝑚 𝑚 𝑛 (𝑝 ) +𝜋5 + ̃𝑚+1 − 𝑝̃𝑚−1 𝑛 + 𝜋6 𝜏 (1.5.2) 𝑛 𝑛 𝑛 𝜌̃𝑔𝑚 2 𝜌̃𝑔𝑚 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 В уравнении (1.5.2) знак «±» обусловлен возможным изменением направления движения газа: в каждой скобке «+» выбирается при положительной компоненте скорости фильтрации газа, умножаемой на данную скобку, а «−» — при отрицательной компоненте [7, 8]. Уравнение импульса в безразмерных переменных 𝜕𝑢̃ 𝜕𝑢̃ 𝜋7 𝜕𝑝̃ 𝑢̃ 𝑇̃𝑔1.5 = −𝜋3 𝑢̃ − − 𝜋8 − 𝜋9 ; 𝜕𝑡̃ 𝜕𝑥̃ 𝜌̃𝑔 𝜕𝑥̃ 𝜌̃𝑔 𝑐𝑠2 + 𝑇̃𝑔 Преобразуем уравнение сохранения импульса в явное конечноразностное уравнение, из которого определим скорость фильтрации газа. Заменяем первые производные по времени и пространству конечными разностями: 20
𝑛+1 𝑛 𝜕𝑢̃ 𝑢̃𝑚 − 𝑢̃𝑚 = ; 𝜕𝑡̃ 𝜏 𝑛 𝑛 𝜕𝑢̃ 𝑢̃𝑚+1 − 𝑢̃𝑚−1 = ; 𝜕𝑥̃ 2ℎ 𝑛 𝑛 𝜕𝑝̃ 𝑝̃𝑚+1 − 𝑝̃𝑚−1 = ; 𝜕𝑥̃ 2ℎ при этом получим алгебраические уравнения: 𝑛 𝑛 𝑛 𝑛+1 𝑛 𝑢̃𝑚 − 𝑢̃𝑚 𝑢̃𝑛 − 𝑢̃𝑚−1 𝜋7 𝑝̃𝑚+1 − 𝑝̃𝑚−1 𝑛 𝑚+1 = −𝜋3 𝑢̃𝑚 − 𝑛 − 𝜋8 𝜏 2ℎ 𝜌̃𝑔𝑚 2ℎ 𝑛 1.5 𝑛 (𝑇 ̃𝑔𝑚 ) 𝑢̃𝑚 − 𝜋9 𝑛 𝑛 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 𝑛+1 Выведем 𝑢̃𝑚 : 𝑛+1 𝑢̃𝑚 = 𝑛 𝑢̃𝑚 𝑛 1.5 𝑛 𝑛 𝑛 𝑛 ) 𝑢̃𝑚+1 − 𝑢̃𝑚−1 𝜋9 𝜏 (𝑇̃𝑔𝑚 𝜋7 𝜏 𝑝̃𝑚+1 − 𝑝̃𝑚−1 − 𝑛 (1 − 𝜋3 𝜏 𝑛 𝑛 )−𝜌 2ℎ 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 ̃𝑔𝑚 2ℎ − 𝜋8 𝜏; 𝑛+1 𝑢̃𝑚 = − 𝑛 𝑢̃𝑚 𝑛 1.5 ) 𝜋3 𝑟ℎ 𝑛 𝜋9 𝜏 (𝑇̃𝑔𝑚 𝑛 (𝑢̃𝑚+1 − 𝑢̃𝑚−1 ) − 𝑛 (1 − 𝑛 ) − 𝜋8 𝜏 2 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 𝜋7 𝑟ℎ 𝑛 𝑛 ). ̃𝑚+1 − 𝑝̃𝑚−1 𝑛 (𝑝 2𝜌̃𝑔𝑚 (1.5.3) Уравнение неразрывности в безразмерных переменных 𝜕 𝑝̃ 𝜕 𝑝̃𝑢̃ ( ) = −𝜋3 ( ); 𝜕𝑡̃ 𝑇̃𝑔 𝜕𝑥̃ 𝑇̃𝑔 Преобразуем уравнение неразрывности в неявное конечно- разностное уравнение, из которого, с учетом уравнения состояния совершенного газа, методом прогонки находится давление газа. Заменяем первые производные по времени и пространству конечными разностями: 21
𝜕 𝑝̃ ( )= 𝜕𝑡̃ 𝑇̃𝑔 𝜕 𝑝̃𝑢̃ ( )= 𝜕𝑥̃ 𝑇̃𝑔 𝑛+1 𝑛 𝑝̃𝑚 𝑝̃𝑚 − 𝑛+1 𝑛 𝑇̃𝑔𝑚 𝑇̃𝑔𝑚 𝜏 ; 𝑛+1 𝑛+1 𝑛+1 𝑛+1 𝑝̃𝑚+1 𝑢̃𝑚+1 𝑝̃𝑚−1 𝑢̃𝑚−1 − 𝑛+1 𝑛+1 𝑇̃𝑔𝑚+1 𝑇̃𝑔𝑚−1 2ℎ ; при этом получим систему алгебраических уравнений 𝑛+1 𝑛 𝑝̃𝑚 𝑝̃𝑚 − 𝑛+1 𝑛 𝑇̃𝑔𝑚 𝑇̃𝑔𝑚 𝜏 𝑛+1 𝑛+1 𝑛+1 𝑛+1 𝑝̃𝑚+1 𝑢̃𝑚+1 𝑝̃𝑚−1 𝑢̃ − ̃ 𝑛+1𝑚−1 𝑛+1 ̃ 𝑇𝑔𝑚+1 𝑇𝑔𝑚−1 = −𝜋3 ; 2ℎ Умножим уравнение на 𝜏 = 𝑟ℎ2 𝑛+1 𝑛+1 𝑛+1 𝑛+1 𝑛+1 𝑛 𝑝̃𝑚 𝑝̃𝑚 𝜋3 𝑟ℎ 𝑝̃𝑚+1 𝑢̃𝑚+1 𝑝̃𝑚−1 𝑢̃𝑚−1 = − ( − ) 𝑛+1 𝑛 𝑛+1 𝑛+1 2 𝑇̃𝑔𝑚 𝑇̃𝑔𝑚 𝑇̃𝑔𝑚+1 𝑇̃𝑔𝑚−1 𝑛+1 𝑝̃𝑚−1 𝑛+1 𝑛+1 𝑛+1 𝑛 𝜋3 𝑟ℎ 𝑢̃𝑚−1 𝑝̃𝑚 𝜋3 𝑟ℎ 𝑢̃𝑚+1 𝑝̃𝑚 𝑛+1 ̃𝑚+1 ( (− 𝑛+1 ) + ̃ 𝑛+1 + 𝑝 𝑛+1 ) = ̃ 𝑛 2 𝑇̃𝑔𝑚−1 2 𝑇̃𝑔𝑚+1 𝑇𝑔𝑚 𝑇𝑔𝑚 𝑛+1 𝑝̃𝑚−1 (− 𝑛+1 𝑛+1 𝑛+1 𝜋3 𝑟ℎ 𝑢̃𝑚−1 𝑝̃𝑚 𝜋3 𝑟ℎ 𝑢̃𝑚+1 𝑛+1 𝑛 , (1.5.4) ) + 𝑛+1 + 𝑝̃𝑚+1 ( ) = 𝜌̃𝑔𝑚 𝑛+1 𝑛+1 2 𝑇̃𝑔𝑚−1 2 𝑇̃𝑔𝑚+1 𝑇̃𝑔𝑚 Каждое разностное уравнение для нахождения давления, содержит 𝑛+1 𝑛+1 𝑛+1 на каждом слое три неизвестных значения 𝑝̃𝑚−1 , 𝑝̃𝑚 , 𝑝̃𝑚+1 , для нахождения которых необходимо решать трёхдиагональную матрицу. Для этого будем использовать метод прогонки. Заменим 𝑐𝑜𝑛𝑠𝑡4 = 𝑛+1 𝑝̃𝑚−1 𝜋3 𝑟ℎ 2 . 𝑛+1 𝑛+1 𝑛+1 𝑢̃𝑚−1 𝑝̃𝑚 𝑢̃𝑚+1 𝑛+1 𝑛 (−𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) + 𝑛+1 + 𝑝̃𝑚+1 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) = 𝜌̃𝑔𝑚 ̃ ̃ ̃ 𝑇𝑔𝑚−1 𝑇𝑔𝑚 𝑇𝑔𝑚+1 При m=2; 𝑝̃1𝑛+1 (−𝑐𝑜𝑛𝑠𝑡4 𝑢̃1𝑛+1 𝑝̃2𝑛+1 𝑢̃3𝑛+1 𝑛+1 ̃3 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) = 𝜌̃2𝑛 𝑛+1 ) + ̃ 𝑛+1 + 𝑝 ̃ 𝑇𝑔1 𝑇𝑔2 𝑇̃𝑔3 22
Перенесем неизвестные влево, известные вправо 𝑝̃2𝑛+1 𝑢̃3𝑛+1 𝑢̃1𝑛+1 𝑛+1 𝑛 𝑛+1 ̃3 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) = 𝜌̃2 − 𝑝̃1 (−𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) 𝑛+1 + 𝑝 𝑇̃𝑔2 𝑇̃𝑔3 𝑇̃𝑔1 𝑏2 𝑝̃2𝑛+1 + 𝑐2 𝑝̃3𝑛+1 = 𝑑2 𝑏2 = 1 𝑛+1 𝑇̃𝑔2 , 𝑢̃3𝑛+1 𝑐2 = (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ), 𝑇̃𝑔3 𝑑2 = 𝜌̃2𝑛 − 𝑝̃1𝑛+1 𝑢̃1𝑛+1 (−𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) 𝑇̃𝑔1 Из (3) выразим 𝑝̃2𝑛+1 : 𝑝̃2𝑛+1 = 𝑐 𝑑2 𝑏2 𝑏2 где α2 = − 2 , β2 = 𝑑2 𝑐2 𝑛+1 − 𝑝̃ = α2 𝑝̃3𝑛+1 + β2 , 𝑏2 𝑏2 3 (1.5.5) . При m=3: 𝑝̃2𝑛+1 (−𝑐𝑜𝑛𝑠𝑡4 𝑢̃2𝑛+1 𝑝̃3𝑛+1 𝑢̃4𝑛+1 𝑛+1 ̃4 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) = 𝜌̃3𝑛 𝑛+1 ) + ̃ 𝑛+1 + 𝑝 ̃ 𝑇𝑔2 𝑇𝑔3 𝑇̃𝑔4 𝑎3 𝑝̃2𝑛+1 + 𝑏3 𝑝̃3𝑛+1 + 𝑐3 𝑝̃4𝑛+1 = 𝑑3 , (1.5.6) где 𝑢̃2𝑛+1 𝑎3 = −𝑐𝑜𝑛𝑠𝑡4 𝑛+1 , 𝑇̃𝑔2 𝑏3 = 1 𝑛+1 𝑇̃𝑔3 , 23
𝑢̃4𝑛+1 𝑐3 = 𝑐𝑜𝑛𝑠𝑡4 𝑛+1 , 𝑇̃𝑔4 𝑑3 = 𝜌̃3𝑛 . Подставляем выражение (1.5.5) в (1.5.6): 𝑎3 (α2 𝑝̃3𝑛+1 + β2 ) + 𝑏3 𝑝̃3𝑛+1 + 𝑐3 𝑝̃4𝑛+1 = 𝑑3 Теперь оно содержит две неизвестных, выразим 𝑝̃3𝑛+1 : 𝑝̃3𝑛+1 = − 𝑐3 𝑑3 − 𝑎3 𝛽2 𝑝̃4𝑛+1 + = 𝛼3 𝑝̃4𝑛+1 + 𝛽3 , 𝑏3 + 𝑎3 𝛼2 𝑏3 + 𝑎3 𝛼2 где 𝛼3 = − 𝛽3 = 𝑐3 , 𝑏3 + 𝑎3 𝛼2 𝑑3 − 𝑎3 𝛽2 . 𝑏3 + 𝑎3 𝛼2 Продолжим вычислять значения 𝛼𝑚 , 𝛽𝑚 пока не дойдем до последнего уравнения. При m=x-3: 𝑛+1 𝑝̃𝑥−4 (−𝑐𝑜𝑛𝑠𝑡4 𝑛+1 𝑛+1 1 𝑢̃𝑥−4 𝑝̃x−3 𝑢̃𝑥−2 𝑛+1 𝑛 ) + 𝑛+1 + 𝑝̃𝑥−2 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) = 𝜌̃𝑥−3 1 𝑇̃𝑔𝑥−4 𝑇̃𝑔𝑥−3 𝑇̃𝑔𝑥−2 Перенесем неизвестные влево, известные вправо 𝑛+1 𝑝̃𝑥−4 (−𝑐𝑜𝑛𝑠𝑡4 𝑛+1 𝑛+1 𝑛+1 𝑢̃𝑥−4 𝑝̃x−3 𝑢̃𝑥−2 𝑛 𝑛+1 ̃𝑥−3 − 𝑝̃𝑥−2 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) 𝑛+1 ) + ̃ 𝑛+1 = 𝜌 𝑇̃𝑔𝑥−4 𝑇𝑔𝑥−3 𝑇̃𝑔𝑥−2 𝑛+1 𝑛+1 𝑎𝑥−3 𝑝̃𝑥−4 + 𝑏𝑥−3 𝑝̃x−3 = 𝑑𝑥−3 где 𝑎𝑥−3 𝑛+1 𝑢̃𝑥−4 = −𝑐𝑜𝑛𝑠𝑡4 𝑛+1 , 𝑇̃𝑔𝑥−4 24
𝑏𝑥−3 = 𝑑𝑥−3 = 𝑛 𝜌̃𝑥−3 − 1 𝑛+1 , 𝑇̃𝑔𝑥−3 𝑛+1 𝑝̃𝑥−2 𝑛+1 𝑢̃𝑥−2 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) 𝑇̃𝑔𝑥−2 𝑛+1 𝑛+1 𝑎𝑥−3 (α𝑥−4 𝑝̃𝑥−3 + β𝑥−4 ) + 𝑏𝑥−3 𝑝̃x−3 = 𝑑𝑥−3 𝑛+1 𝑝̃x−3 = βx−3 = 𝑑𝑥−3 − 𝑎𝑥−3 β𝑥−4 𝑏𝑥−3 + 𝑎𝑥−3 α𝑥−4 𝑛+1 𝑛+1 𝑛 𝑛+1 Зная 𝑝̃x−3 , можно найти 𝑝̃𝑥−4 = 𝛼𝑥−4 𝑝̃𝑥−3 + 𝛽𝑥−4 . Зная 𝑝̃x−4 , можно 𝑛+1 𝑛 найти 𝑝̃𝑥−5 = 𝛼𝑥−5 𝑝̃𝑥−4 + 𝛽𝑥−5 . С каждым шагом узнаем значение новой переменной, номер которой на 1 меньше предыдущей. Так добираемся до m = 2. Все переменные найдены, задача решена. Плотность газа рассчитывается из уравнения состояния совершенного газа. Таким образом, система уравнений (1.4.19) преобразуется: 𝑛 𝑛 ); 𝑛 𝑇𝑚𝑛+1 = 𝑇𝑚𝑛 (1 − 𝜋1 𝜏) + 𝜋1 𝜏𝑇𝑔𝑚 + 𝜋2 𝑟(𝑇𝑚+1 − 2𝑇𝑚𝑛 + 𝑇𝑚−1 𝑛+1 𝑛 𝑇̃𝑔𝑚 = 𝑇̃𝑔𝑚 (1 − 𝜋4 𝜏 𝜋3 𝑟ℎ 𝑛 𝜋4 𝜏 𝑛 𝑛 𝑛 𝑛 ̃ − 𝑢 ̃ ∓ 4𝑇 ± 𝑇 + ) (±3𝑇 ) 𝑚 𝑔𝑚 𝑔𝑚∓1 𝑔𝑚∓2 𝑛 𝑛 𝑇𝑐𝑚 + 𝜌̃𝑔𝑚 2 𝜌̃𝑔𝑚 1.5 𝑛 𝑛 𝑛−1 𝑛 𝑛 )2 (𝑇̃𝑔𝑚 ) (𝑢̃𝑚 𝑝̃𝑚 − 𝑝̃𝑚 𝜋5 𝜋3 𝑟ℎ 𝑛 𝑢̃𝑚 𝑛 (𝑝̃𝑚+1 − 𝑝̃𝑚−1 ) 𝑛 + 𝜋6 𝜏 +𝜋5 + 𝑛 𝑛 𝑛 𝜌̃𝑔𝑚 2 𝜌̃𝑔𝑚 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 𝑛+1 𝑢̃𝑚 = − 𝑛 𝑢̃𝑚 𝑛 1.5 ) 𝜋3 𝑟ℎ 𝑛 𝜋9 𝜏 (𝑇̃𝑔𝑚 𝑛 (𝑢̃𝑚+1 − 𝑢̃𝑚−1 ) − 𝑛 (1 − 𝑛 ) − 𝜋8 𝜏 2 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 𝜋7 𝑟ℎ 𝑛 𝑛 ); ̃𝑚+1 − 𝑝̃𝑚−1 𝑛 (𝑝 2𝜌̃𝑔𝑚 𝑛+1 𝑝̃𝑚−1 (− (1.5.7) 𝑛+1 𝑛+1 𝑛+1 𝜋3 𝑟ℎ 𝑢̃𝑚−1 𝑝̃𝑚 𝜋3 𝑟ℎ 𝑢̃𝑚+1 𝑛+1 𝑛 ; ) + 𝑛+1 + 𝑝̃𝑚+1 ( ) = 𝜌̃𝑔𝑚 𝑛+1 𝑛+1 2 𝑇̃𝑔𝑚−1 2 𝑇̃𝑔𝑚+1 𝑇̃𝑔𝑚 25
𝑛+1 𝜌𝑚 𝑛+1 𝑝𝑚 = 𝑛+1 ; 𝑇𝑔𝑚 Для удобства описания алгоритмов произведения констант запишем в новые постоянные величины: 𝑐𝑜𝑛𝑠𝑡1 = −𝜋1 𝜏, 𝑐𝑜𝑛𝑠𝑡11 = 1 + 𝑐𝑜𝑛𝑠𝑡1, 𝑐𝑜𝑛𝑠𝑡2 = 𝜋2 𝑟, 𝑐𝑜𝑛𝑠𝑡3 = 𝜋4 𝜏, 𝑐𝑜𝑛𝑠𝑡4 = 𝜋3 𝑟ℎ , 2 𝑐𝑜𝑛𝑠𝑡5 = 𝜋6 𝜏, 𝑐𝑜𝑛𝑠𝑡6 = 𝜋9 𝜏, 𝑐𝑜𝑛𝑠𝑡7 = 𝜋8 𝜏, 𝑐𝑜𝑛𝑠𝑡8 = 𝜋7 𝑟ℎ , 2 𝑐𝑜𝑛𝑠𝑡9 = 𝜋5 𝑐𝑜𝑛𝑠𝑡4, Тогда система уравнений (1.5.5) примет вид: 𝑛 𝑛 ) 𝑛 𝑇𝑚𝑛+1 = 𝑇𝑚𝑛 𝑐𝑜𝑛𝑠𝑡11 − 𝑐𝑜𝑛𝑠𝑡1 𝑇𝑔𝑚 + 𝑐𝑜𝑛𝑠𝑡2(𝑇𝑚+1 − 2𝑇𝑚𝑛 + 𝑇𝑚−1 ; 𝑛+1 𝑛 𝑇̃𝑔𝑚 = 𝑇̃𝑔𝑚 (1 − 𝑐𝑜𝑛𝑠𝑡3 𝑛 𝑛 𝑛 𝑛 ∓ 4𝑇𝑔𝑚∓1 ± 𝑇𝑔𝑚∓2 ) − 𝑐𝑜𝑛𝑠𝑡4 𝑢̃𝑚 (±3𝑇𝑔𝑚 ) 𝑛 𝜌̃𝑔𝑚 𝑛 𝑛−1 𝑛 𝑐𝑜𝑛𝑠𝑡3 𝑛 𝑝̃𝑚 − 𝑝̃𝑚 𝑢̃𝑚 𝑛 𝑛 ̃ + 𝑛 𝑇𝑐𝑚 + 𝜋5 + 𝑐𝑜𝑛𝑠𝑡9(𝑝̃𝑚+1 − 𝑝̃𝑚−1 ) 𝑛 𝑛 𝜌̃𝑔𝑚 𝜌̃𝑔𝑚 𝜌̃𝑔𝑚 𝑛 1.5 𝑛 )2 (𝑇̃𝑔𝑚 ) (𝑢̃𝑚 + 𝑐𝑜𝑛𝑠𝑡5 ; 𝑛 𝑛 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 26
𝑛+1 𝑢̃𝑚 𝑛 1.5 ̃𝑔𝑚 (𝑇 ) 𝑐𝑜𝑛𝑠𝑡6 𝑛 𝑛 𝑛 )− 𝑛 = 𝑢̃𝑚 − 𝑢̃𝑚−1 (1 − 𝑐𝑜𝑛𝑠𝑡4(𝑢̃𝑚+1 𝑛 ) − 𝑐𝑜𝑛𝑠𝑡7 𝜌̃𝑔𝑚 𝑐̃𝑠2 + 𝑇̃𝑔𝑚 − 𝑛+1 𝑝̃𝑚−1 𝑐𝑜𝑛𝑠𝑡8 𝑛 𝑛 (𝑝̃𝑚+1 − 𝑝̃𝑚−1 ) − 𝜔1 ; 𝑛 𝜌̃𝑔𝑚 (1.5.6) 𝑛+1 𝑛+1 𝑛+1 𝑢̃𝑚−1 𝑝̃𝑚 𝑢̃𝑚+1 𝑛+1 𝑛 − 𝜔2 ; (−𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) + 𝑛+1 + 𝑝̃𝑚+1 (𝑐𝑜𝑛𝑠𝑡4 𝑛+1 ) = 𝜌̃𝑔𝑚 𝑇̃𝑔𝑚−1 𝑇̃𝑔𝑚 𝑇̃𝑔𝑚+1 𝑛+1 𝜌𝑚 𝑛+1 𝑝𝑚 = 𝑛+1 ; 𝑇𝑔𝑚 где 𝜔1 = 1 𝑛 𝑛 𝑛 𝑛 𝑛 ). ̃ 𝑚+2 − 4𝑢̃𝑚+1 + 6𝑢̃𝑚 − 4𝑢̃𝑚−1 + 𝑢̃𝑚−2 ∗ (𝑢 𝜔1 𝜔2 = 1 𝑛 𝑛 𝑛 𝑛 𝑛 (𝑝̃ ). − 4𝑝̃𝑚+1 + 6𝑝̃𝑚 − 4𝑝̃𝑚−1 + 𝑝̃𝑚−2 𝜔2∗ 𝑚+2 Здесь 𝜔1 , 𝜔2 − демпфирующие члены, необходимые для уменьшения дисперсионной ошибки. Такая ошибка может возникать из-за нелинейности системы уравнений, которая в данном случае может приводить к колебаниям получаемого решения, при этом сами демпфирующие члены не изменяют формальную точность решения [7,8]. Алгоритм решения системы уравнений. Конечными-разностями решаем первое – третье уравнения системы (1.5.6), далее методом прогонки решаем четвертое уравнение и затем тривиально решаем пятое. Также необходимо задать значения давления и температуры газа в фиктивных точках, расположенных за пределами пористого элемента. 27
Глава 2. Результаты решения 2.1. Тестовая задача Протестируем правильность работы алгоритмов, решив одномерную задачу о нестационарном режиме охлаждения газом, протека0ющим через твёрдую пористую неподвижную среду, в которой происходит процесс саморазогрева пористого элемента. Подобное решение задачи было предложено для описания процесса охлаждения аварийного реактора на четвертом энергоблоке ЧАЭС [6-8,10]. Дополним уравнение энергии твердой фазы слагаемым, отвечающим за тепловыделение в твердой фазе 𝑄0 (1 − 𝑎)𝐶: 𝜕𝑇𝑐 𝜕 2 𝑇𝑐 (1 − 𝑎)𝜌𝑐 𝐶𝑐 = −𝛼(𝑇𝑐 − 𝑇𝑔 ) + 𝑄0 (1 − 𝑎)𝐶 + (1 − 𝑎)𝜆𝑐 , (2.1.1) 𝜕𝑡 𝜕𝑥 2 где 𝑄0 – константа, определяющая интенсивность тепловыделения, Дж ; м3 ∙с 𝐶 – относительная концентрация реагирующего вещества; Скорость убывания концентрации реагирующего вещества прямо пропорциональна самой концентрации: 𝜕𝐶 = 𝑘2 𝐶, 𝜕𝑡 (2.1.2) где 𝑘2 – коэффициент, определяющий уменьшение тепловыделения, с−1 ; Расписав в полную производную (2.1.2) и решив при условии, что в начальный момент времени С|𝑡=0 = 1, получим выражение для концентрации реагирующего вещества: С = 𝑒𝑥𝑝(−𝑘2 𝑡). (2.1.3) 28
Подставив (2.1.2) в (2.1.1) и обезразмерив полученное уравнение выведем уравнение притока тепла в безразмерных переменных: 𝑛 𝑇𝑚𝑛+1 = 𝑇𝑚𝑛 (1 − 𝜋1 𝜏) + 𝜋1 𝜏𝑇𝑔𝑚 + 𝑄𝑒𝑥𝑝(−𝜀𝜏) 𝑛 𝑛 ), + 𝜋2 𝑟(𝑇𝑚+1 − 2𝑇𝑚𝑛 + 𝑇𝑚−1 (2.1.4) с параметрами подобия: 𝑄= 𝑄0 𝑡∗ , 𝜌𝑐 𝐶𝑐 𝑇∗ 𝜀 = 𝑘2 𝑡∗ Решения системы (1.4.6) с учетом преобразований (2.1.1-2.1.4) и краевыми условиями получены из (1.3.5) при следующих значениях: 𝑝0 = 1.5; 𝑇𝑔0 = 1; 𝑝ℎ = 1; (2.1.5) Для рассматриваемого решения системы безразмерные параметры примут значения: 𝜋1 = 7.06 ∙ 10−4 ; 𝜋2 = 5.93 ∙ 10−9 ; 𝜋3 = 0.333; 𝜋4 = 2.778; 𝜋5 = 2.34 ∙ 10−2 ; 𝜋6 = 1.85 ∙ 103 ; 𝜋7 = 2.179; 𝜋8 = 4.68 ∙ 102 , 𝑄 = 1.65 ∙ 10−4 ; 𝜀 = 10−7 ; Они соответствуют следующим размерным величинам: 𝐻 = 10 м; 𝑡∗ = 1 𝑐; 𝑢∗ = 1 𝑝∗ = 105 Па; 𝜌∗ = 1.2 м ; 𝑇 = 300 К; с ∗ кг кг Дж 3 2 ; 𝜌 = 2.2 ∙ 10 ; 𝐶 = 9.2 ∙ 10 ; м3 𝑐 м3 𝑐 кг · К 𝛼 = 103 Дж Дж м 3 ; 𝐶 = 10 ; 𝑔 = 9.8 ; 𝑘1 = 10−8 м2 ; 𝑔𝑝 3 2 м ·К∙с кг · К с 𝑎∗ = 10 Дж ; с = 1.458 ∙ 10−8 м2 ; с𝑠2 = 110.4 К; м2 ∙ К ∙ с 𝑠1 (2.1.3) 29
𝑅 = 8.314 Дж Вт ; 𝜆 = 1.2 ; 𝑎 = 0.3; 𝑘2 = 10−7 с−1 ; (кг · К) м·К 𝑄0 = 105 Дж ; м3 ∙ с На рис 2.1.1-2.1.5 представлено решение системы уравнений (2.1.1) с краевыми условиями (2.1.2) и с набором значений размерных величин (2.1.3). Более подробно о решении и результатах, представленных на графиках, описано в статье [7]. Рис. 2.1.1 Температура твердой фазы 30
Рис. 2.1.2 Температура газа Рис. 2.1.3 Скорость фильтрации газа 31
Рис. 2.1.4 Давление газа Рис. 2.1.5 Плотность газа 32
Сравнение полученных результатов тестовой задачи с результатами, полученными в [7], показывает, что реализованный алгоритм для решения математической модели, представленной в виде системы уравнений (1.5.6) с краевыми условиями (1.4.4), работает корректно. 2.2. Результаты вычислительных экспериментов о нестационарном режиме нагревания пористого элемента Для определения основных закономерностей течения горячего газа через пористый объект проанализируем решение системы (1.4.6) с краевыми условиями (1.4.4). В качестве начальных условий выбиралось стационарное решение из (1.3.1). Пусть давление на входе в пористый объект резко достигает 1.5 атм и после этого остается постоянным. На выходе из пористого элемента газ попадает в свободное пространство с постоянным давлением 1 атм. Температура газа на входе составляет 400 К и не меняется. 𝑝0 = 1.5 атм; 𝑇𝑔0 = 400 К; 𝑝ℎ = 1 атм; (2.2.1) Пусть свойства твердой среды близки к свойствам бетона, а параметры газовой фазы соответствуют воздуху. Тогда для рассматриваемого решения системы безразмерные параметры примут значения: 𝜋1 = 7.06 ∙ 10−4 ; 𝜋2 = 5.93 ∙ 10−9 ; 𝜋3 = 0.333; 𝜋4 = 2.778; 𝜋5 = 2.34 ∙ 10−2 ; 𝜋6 = 1.85 ∙ 103 ; 𝜋7 = 2.179; 𝜋8 = 4.68 ∙ 102 . Они соответствуют следующим размерным величинам: 𝐻 = 10 м; 𝑡∗ = 1 𝑐; 𝑢∗ = 1 𝑝∗ = 105 Па; 𝜌∗ = 1.2 м ; 𝑇 = 300 К; с ∗ кг кг Дж 3 2 ; 𝜌 = 2.2 ∙ 10 ; 𝐶 = 9.2 ∙ 10 ; м3 𝑐 м3 𝑐 кг · К 33
𝛼 = 103 Дж Дж м 3 ; 𝐶 = 10 ; 𝑔 = 9.8 ; 𝑘 = 10−8 м2 ; м3 · К ∙ с 𝑔𝑝 кг · К с2 1 𝑎∗ = 10 Дж ; с𝑠1 = 1.458 ∙ 10−8 м2 ; с𝑠2 = 110.4 К; 2 м ∙К∙с 𝑅 = 8.314 (2.2.2) Дж Вт ; 𝜆 = 1.2 ; 𝑎 = 0.3. (кг · К) м·К На рисунках 2.2.1-2.2.5 показано распределение температуры фаз, плотности газа, давления, скорости фильтрации по высоте пористого объекта. Через четыре часа после начала течения горячего газа твердая фаза достаточно прогрелась. Однако для перехода системы к стационарному режиму, когда температура твердой среды стала равной температуре газа, потребовалось чуть более пяти часов. Рис. 2.2.1 Температура твердой фазы 34
Из рисунка 2.2.1 видно, что температура твердой фазы на входе в пористый объект быстро возрастет, что приводит к быстрому изменению всех остальных параметров. Рост температуры газа со временем приводит к уменьшению плотности газа. При разогреве, газу становится сложнее проходить через пористый объект, и скорость фильтрации уменьшается. В сочетании с уменьшением перепада температур это приводит к тому, что с течением времени скорость изменения температуры уменьшается. Рис. 2.2.2 Температура газа 35
Рис. 2.2.3 Плотность газа На рисунке 2.2.3 кривая для 10 минут показывает, как в первое время происходит возрастание плотности газа (достигая максимального значения 1.34 кг/м3 при x = 2.5 м), а затем его уменьшение. Данный эффект связан с интенсивным уменьшение температуры газа (кривая температуры газа на начальном этапе имеет выпуклость вниз). Газ, доходя до x = 2.5 м, отдает большую часть своей тепловой энергии пористой среде, при этом давление уменьшается лишь на 25%, что и приводит к возрастанию плотности. Это явление сохраняется до тех пор, пока изменение температуры газа вдоль пористого элемента не придет к более плавному виду. 36
Рис. 2.2.4 Давление Рис. 2.2.5 Скорость фильтрации газа 37
На рисунке 2.2.5 показано распределение скорости фильтрации газа, из которого видно, что кривая для t = 10 минут сначала убывает до значения 2.34 м/с, а затем идет на увеличение. При этом минимум скорости находится в той же точке пористого элемента, что и максимум плотности газа. С течением времени средняя скорость фильтрации газа уменьшается, однако при рассмотрении кривой для t = 4 часов видно, что скорость фильтрации вдоль пористого элемента возрастает быстрее, чем при t = 20 мин или 1 часу. Всё это указывает на то, что плотность газа существенно влияет на скорость фильтрации. На рисунке 2.2.6 показано изменение со временем температур твердой и газовой фаз на выходе из пористого элемента. Рис. 2.2.6 Изменение температур фаз на выходе 38
Рассмотрим влияние пористости на нестационарный режим нагрева пористого объекта. На рисунках 2.2.7 и 2.2.9 представлены графики изменения температуры твердой и газовой фазы для различной пористости при одинаковом моменте времени, а на рисунках 2.2.8 и 2.2.10 – разница температур фаз на выходе из пористого элемента при различной пористости. Здесь и в дальнейшем под разностью температур понимается выражение 𝑇|𝑎𝑖+1 − 𝑇|𝑎𝑖 , при t = 30 минут, где 𝑎𝑖 – пористость 𝑎𝑖 ∈ (0.9: 0.1) . Рис. 2.2.7 Температура твердой фазы 39
. Рис. 2.2.8 Разница температур твердой фазы на выходе из пористого элемента при различной пористости Рис. 2.2.9 Температура газа 40
Рис. 2.2.10 Разница температур газа на выходе из пористого элемента при различной пористости Пористость напрямую связана с площадью поверхности твердой фазы. Из увеличения пористости следует увеличение площади соприкосновения твердой фазы с газовой. Также с увеличением пористости растет расход протекающего через элемент газа и уменьшается масса твердой среды в пористом элементе. За счет этого общая скорость теплообмена и нагрев твердой фазы увеличиваются. 2.3. Результаты вычислительных экспериментов о нестационарном режиме охлаждения пористого элемента Рассмотрим результаты вычислений для течения холодного газа через неподвижный разогретый пористый объект. Пусть давление на входе в пористый объект резко достигает 1.5 атм и после этого остается постоянным. На выходе из пористого элемента газ попадает в свободное 41
пространство с постоянным давлением 1 атм. Температура газа на входе составляет 300 К и не меняется. 𝑝0 = 1.5 атм; 𝑇𝑔0 = 300 К; 𝑝ℎ = 1 атм; (2.3.1) Пусть свойства твердой среды близки к свойствам бетона, а параметры газовой фазы соответствуют воздуху. Тогда для рассматриваемого решения системы безразмерные параметры примут значения: 𝜋1 = 7.06 ∙ 10−4 ; 𝜋2 = 5.93 ∙ 10−9 ; 𝜋3 = 0.333; 𝜋4 = 2.778; 𝜋5 = 2.34 ∙ 10−2 ; 𝜋6 = 1.85 ∙ 103 ; 𝜋7 = 2.179; 𝜋8 = 4.68 ∙ 102 . Они соответствуют следующим размерным величинам: 𝐻 = 10 м; 𝑡∗ = 1 𝑐; 𝑢∗ = 1 𝑝∗ = 105 Па; 𝜌∗ = 1.2 м ; 𝑇 = 300 К; с ∗ кг кг Дж 3 2 ; 𝜌 = 2.2 ∙ 10 ; 𝐶 = 9.2 ∙ 10 ; 𝑐 𝑐 м3 м3 кг · К 𝛼 = 103 Дж Дж м 3 ; 𝐶 = 10 ; 𝑔 = 9.8 ; 𝑘 = 10−8 м2 ; м3 · К ∙ с 𝑔𝑝 кг · К с2 1 𝑎∗ = 10 Дж ; с𝑠1 = 1.458 ∙ 10−8 м2 ; с𝑠2 = 110.4 К; 2 м ∙К∙с 𝑅 = 8.314 (2.3.2) Дж Вт ; 𝜆 = 1.2 ; 𝑎 = 0.3. (кг · К) м·К На рисунках 2.3.1-2.3.5 показано распределение температуры фаз, плотности газа, давления, скорости фильтрации по высоте пористого объекта. Через пять часов после начала течения холодного газа, пористый объект полностью остыл, и система пришла к стационарному режиму. 42
Рис. 2.3.1 Температура твердой фазы Из рисунков 2.3.1 – 2.3.2 видно, что температура твердой фазы на входе в пористый объект резко убывает. Как и было описано ранее, наибольшая скорость нагрева (остывания) газа пористым элементом имеется в первое время, пока разница между температурами фаз достаточно велика. Со временем разность температур уменьшается, вместе с этим уменьшается скорость нагрева, но увеличивается время поддержания температурного диапазона низких температур. 43
Рис. 2.3.2 Температура газа Ближе к входу в пористый элемент, скорость фильтрации меньше, чем на выходе из него. В совокупности с тем, что температура твердой фазы еще достаточно велика, это приводит к тому, что кривые распределения температуры газа в промежутке от 0 секунд до 40 минут имеют выпуклость вверх, после – выпуклость вниз. 44
Рис. 2.3.3 Плотность газа При режиме нагрева холодного газа пористым элементом возрастание, а затем убывание (или наоборот) распределение плотности газа, как на рисунке 2.2.3, не проявилось. Со сменой выпуклости кривой плотности газ становится менее разряженным, так как из-за постоянного оттока тепла твердая фаза нагревает газ менее эффективно, вплоть до полного остывания системы (𝑇𝑔 = 𝑐𝑜𝑛𝑠𝑡), тогда изменение плотности газа вдоль пористого объекта зависит исключительно от перепада давления. 45
Рис. 2.3.4 Давление Рис. 2.3.5 Скорость фильтрации газа 46
На рисунке 2.3.6 показано изменение со временем температур твердой и газовой фаз на выходе из пористого элемента. Рис. 2.3.6 Изменение температур фаз на выходе Рассмотрим влияние пористости на нестационарный режим охлаждения пористого объекта. На рисунках 2.2.7 и 2.2.9 представлены графики изменения температуры твердой и газовой фазы для различной пористости после 30 минут от начала процесса, а на рисунках 2.3.8 и 2.2.10 – разница температур фаз на выходе из пористого элемента при различной пористости. 47
Рис. 2.3.7 Температура твердой фазы Средняя разница между кривыми температур твердой фазы на входе составляет 4.6 К. При этом на выходе из пористого элемента каждая последующая разница температур, начиная с a=0.1, в среднем в полтора раза меньше предыдущей. 48
Рис. 2.3.8 Разница температур твердой фазы на выходе из пористого элемента при различной пористости Рис. 2.3.9 Температура газа 49
Подобная тенденция видна и с температурами газа (рис 2.2.8 и 2.2.10). График зависимости разницы температур от пористости имеет экспоненциальный характер. Рис. 2.3.10 Разница температур газа на выходе из пористого элемента при различной пористости 50
Заключение Предложена модель и численный метод для моделирования движения газа через пористый теплоаккумулирующий объект. Разработано программное средство, реализующее предложенный численный метод, и проведено его тестирование. Получено численное решение системы уравнений, описывающей процессы движения газа через пористый теплоаккумулирующий объект. Исследованы два случая одномерного нестационарного движения газа при заданном перепаде давления на открытых границах объекта и заданной температуре входящего газа: 1) движение горячего газа через слой холодного теплоаккумулирующего пористого объекта; 2) движение холодного газа через слой разогретого теплоаккумулирующего пористого объекта. Проанализировано влияние пористости твердой среды на распределение температур. Выявлено, что для более энергоемкого теплоаккумулятора необходимо использовать маленькую пористость. В случае необходимости его быстрой зарядки во вред энергоемкости требуется выбрать большую пористость. 51
Литература 1. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987. 636 с. 2. Ентов В.М., Теория фильтрации// Соросовский образовательный журнал. 1998. № 2. с. 121-123. 3. Жуковский Н.Е. О влиянии давления на насыщенные водою пески. // Собр. соч., т. 7. – М.: Гостехиздат, 1950. – с. 73-89. 4. Левин В.А., Луценко Н.А., Фецов С.С. Моделирование движения газа через слой гранулированного теплоаккумулирующего материала с фазовым переходом // Доклады Академии наук. 2018. Т. 479, № 4. C. 386-389 5. Лойцянский Л.Г. Механика жидкости и газа. – Изд. 6-е, перераб. и доп. – М.: Наука, 1987. – 840 с. 6. Луценко Н.А. Одномерный стационарный режим фильтрации газа через слой неподвижного тепловыделяющего конденсированного материала// Дальневосточный мат. журнал. 2002. Т. 3. № 1. с. 123130. 7. Луценко, Н. А. Нестационарные режимы охлаждения пористого тепловыделяющего элемента / Н. А. Луценко // Математическое моделирование. – 2005. – Т. 17. – № 3. – С. 120–128. 8. Луценко, Н. нестационарных А. Численное течений газа моделирование через пористые трехмерных объекты с источниками энерговыделения / Н. А. Луценко // Вычислительная механика сплошных сред. – 2016. – Т. 9. – № 3. – С. 331–344. 9. Маскет М. Течение однородных жидкостей в пористой среде. – М. – Л.: Гостоптехиздат, 1949. 628 с. 10.Маслов В.П., Мясников В.П., Данилов В.Г. Математическое моделирование аварийного блока Чернобыльской АЭС. М.: Наука, 1987. 120 с. 52
11. Нигматулин Р.И. Основы механики гетерогенных сред. М.: Наука, 1978. 336 с. 12. Нигматулин Р.И., Динамика многофазных сред. Т. 1.: Наука, 1987. 464 с. 13. Ольховский, Г. Г. Воздушно-аккумулирующие газотурбинные электростанции (ВАГТЭ) / Г. Г. Ольховский и соавт. – Ижевск: ИКИ, 2011. – 360 c. 14. Павловский Н.Н. О фильтрации воды через земляные плотины. – Л.: Кубуч, 1931. – 259 с. 15. Полубаринова-Кочина П.Я. Теория движения грунтовых вод. М.: Наука, 1977. 664 с. 16. Седов Л.И. Механика сплошной среды. Т.1. М.: Наука, 1976. 536 с. 17. Эглит М.Э. Лекции по основам механики сплошных сред. – М.: Издво Моск. ун-та, 2008. – 318 с. 18. Barbour E. Adiabatic compressed air energy storage with packed bed thermal energy storage. Applied Energy. 2015. Vol. 155. Pp. 804-815. 19. Darcy H. Les fontaines publiques de la ville de Dijon. – Paris, 1856. – 647 p. 20. Luo, X. Overview of current development in compressed air energy storage technology. Energy Procedia. 2014. Vol. 62. Pp. 603-611. 21. Peng H. Modeling on heat storage performance of compressed air in a packed bed system. Applied Energy. 2015. Vol. 160. Pp. 1-9. 53
Отзывы:
Авторизуйтесь, чтобы оставить отзыв