Mathematical modeling of unsteady heat processes in countercurrent heat exchangers.pdf В энергетике, металлургии, химической промышленности и других отраслях широко используются теплообменные аппараты, обычно - противоточные. От качества работы таких аппаратов существенно зависит качество работы других объектов потребителей теплоносителей и всей промышленности в целом. Для проектирования управления теплообменными аппаратами в переходных режимах необходимо знание достаточно точных математических моделей нестационарных теплофизических процессов, протекающих в них. Кроме получения достоверных математических моделей нестационарных противоточных процессов возникает и задача их совместного численного решения с учетом реального изменения теплофизических параметров по длине потоков. 1. Состояние вопроса Переходные теплофизические процессы описываются дифференциальными уравнениями в частных производных. К сожалению, на практике от прямого решения таких уравнений для моделирования переходных процессов в противоточных теплообменных аппаратах отказываются - вводят существенные ограничения на состояние процессов и переходят от дифференциальных уравнений к передаточным функциям [1-3]. Такая модель хоть и описывает температурное поле внутри тепло-обменного аппарата, но формулируется с существенными допущениями, предполагающими постоянство плотности скорости потока и других параметров теплоносителей по длине теплообменника. Рассмотрим математическую модель одномерного нестационарного теплового потока, состоящую из уравнений теплопереноса, движения, неразрывности и состояния [4, 5]: ST ST X d2T q Su Su 1 Sp n S2 и - + u----- = --, - + u- =---- +--- St Sx cp p Sx CpP St дх p Sx p Sx (1) g+ ^ =0, ®(p,p,T)=0. St Sx Здесь T(x,т) - температура теплоносителя, t e[0,T] - время, x ха, хь J - пространственная координата, u( x, t) - скорость потока, p( x, t) - плотность теплоносителя, p( x, t) - давление в потоке, X -коэффициент теплопроводности, cp - удельная теплоемкость при постоянном давлении, q = q + q от внешних источников q и возможных фазовых переходов q . Вт/ удельный тепловой поток м П - коэффициент динамической вязкости, Ф - уравнение состояния неизотермических процессов. Скорости хладагентов в промышленных теплообменных аппаратах достигают десятков метров в секунду. Конвективный теплообмен здесь существенно преобладает над диффузионным, поэтому диффузионным членом в первом уравнении пренебрегают. Во втором уравнении присутствует член, характеризующий вязкость потока. Он определяет диссипативные процессы, приводящие к потере гидродинамического напора по длине теплообменно-го аппарата и, соответственно, - к некоторому нагреванию хладагентов. В реальных аппаратах потеря напора в теплообменниках составляет примерно 1%. Это - незначительная величина, поэтому далее вязкостью и изменением гидродинамического напора, связанного с диссипацией, мы будем пренебрегать. Естественно, будем пренебрегать и теплотой образующейся от диссипации энергии. Именно поэтому в систему (1) мы заранее не добавляли уравнение энергии. Второе и третье уравнения представляют собой систему гиперболических уравнений, которая описывает волновые гидродинамические процессы в сжимаемых средах, когда изменение давления p(x, т) приводит к изменению плотности р(х,т). Как известно, скорость волн определяется в основном скоростью звука в среде, которая в промышленных хладагентах на несколько порядков выше скорости движения среды. В то же время волны повышения и понижения давления переносят импульс, а не сами частицы среды, как это происходит при конвективном теплообмене. Таким образом, в среднем по длине потока влиянием быстротечных волновых процессов повышения и понижения давления на относительно медленный тепловой конвективный перенос можно пренебречь. В итоге мы можем считать, что среда с хладагентом несжимаема и давление в среде по всей длине хладагента постоянно, т.е. p(х,т) = р(т). Невязкость и несжимаемость сред приводят к тому, что второе уравнение в системе (1) исчезает, а вся система с учетом предыдущих допущений приобретает вид: (2) f+Uf=-i- [q+q], f^ 0, p (T)p = 0. ОТ ox cpp ox lp Последнее уравнение состояния Ф вырождается в зависимость плотности среды от ее температуры при заданном давлении [6]. 2. Математическое моделирование противоточных сред На рис. 1 показана принципиальная схема работы противоточного теплообменного аппарата, в котором прямым потоком может быть, например, воздух, а обратным - жидкий кислород. Будем отмечать индексом 1 прямоточную среду (воздух) и индексом 2 - противоточную среду (кислород). Рис. 1. Схема противоточного теплообменного аппарата Fig. 1. Diagram of a counterflow heat exchanger Одномерная постановка задачи предполагает достаточно интенсивное турбулентное перемешивание, такое что температура и скорость в поперечном сечении каждого потока элементарной толщины dx постоянны. Тогда процесс передачи тепла между участками длиной dx через площадь их взаимодействия dF = Hdx (H - ширина разделительной стенки потоков), согласно законам теплопроводности и теплоотдачи, можно описать уравнениями мощностей элементарных тепловых потоков dP [7, 8]: dPj =-ajdF(Ti -T), dPi2 = --Ц2dF(T -T2), dP2 =-a2dF(f2 -T2), где dP1 - мощность теплового потока от теплоносителя 1 к стенке, dP12 - мощность потока в стенке, dP2 - мощность теплового потока от стенки к теплоносителю 2, а1, а2 - коэффициенты теплоотдачи от теплоносителя 1 к стенке и от стенки к теплоносителю 2 соответственно, Т12 - коэффициент теплопроводности стенки, 5 - толщина стенки, T, T2 - температуры стенки со стороны 1-го и 2-го теплоносителей. Здесь T1 > T1 > T2 > T2 . Мы видим, что тепло из среды 1 в среду 2 попадает с некоторым термическим сопротивлением 1 5 1 Rk = ~ + V + а1 Т12 а2 Поскольку в теплообменных аппаратах разделительную стенку делают малой толщины и из высоко теплопроводящих материалов, то значение ее термического сопротивления 5/Х12 « 0, и им пренебрегают. Кроме того, подразумевая качественную внешнюю теплоизоляцию теплообменных аппаратов, тепловыми потерями в окружающую среду также обычно пренебрегают. Тогда теплообмен для каждой среды можно представить как удельную мощность (на единицу объема) теплового потока Вт/ м q = ±kdF (T - T2 ) / dV, где k = (1/а + 1/а2) 1 - коэффициент теплопередачи Вт/ dV = Sdx - объем элементарной м2К ячейки, принимающей (знак плюс) или отдающей (знак минус) тепло, S - площадь поперечного сечения потока. Если течения сопровождаются фазовыми переходами, то необходимо учитывать соответствующие выделения и поглощения тепла. Например, охладитель в виде жидкого кислорода внутри теплообменника всегда закипает и выходит в виде газа. При парообразовании кислород поглощает удельный тепловой поток q2 = q1 из первой противоточной среды (воздух) и удерживает температуру на уровне кипения T2 = T2rm вплоть до полного испарения. Здесь f-p2L2 / Дг2, если есть кипение, 0, если нет кипения, 12 =q1 =-kdF (T1 - T2) / dV, где L2 - удельная теплота парообразования кислорода рода до полного испарения. Таким образом система (2) принимает вид: Дж/ , Дт2 - время кипения ячейки кислокг -kH(T1 - T2 ) + q1 S1 h 1 - + U дх dT^ дх ср1р1 д?2 дт 1 kS~ (T1 - T2 ) + q2 (3) --Un ср 2Р2 dp1U1 = 0, = 0, Р1 (T1 )| P = 0, Р2 (T2 )| P = 0. дх 2 дх Если для какого-либо хладагента отсутствуют табличные данные зависимости р (T) в диапазоне рабочих давлений теплообменника, то необходимо воспользоваться уравнением состояния неидеальных сред (при высоких давлениях p), например в виде уравнения Битти-Бриджмена [9]. Запишем его для одного моля вещества (1 - s)(V + 5) A RT------ - p = 0, (4) V2 V где R - универсальная газовая постоянная, V (T) - объем моля вещества массой M, коэффициенты A = A (l - aV), В = B0 (l - by-), s = С/Vtз, A0, B0, a, b, c - эмпирические постоянные [10]. Объем V из нелинейного уравнения (4) при заданном давленииp можно найти итерационно. Обычно это делается методом Ньютона. После чего мы получаем плотность р (T ) = M/v (T). (5) V (T )■ Для вычисления удельной теплоемкости ср (T) обеих сред используют табличные данные, например [11]. На основе этих данных делается интерполяция функции ср (T) . В системе (3) первые два уравнения конвективного теплопереноса нуждаются в граничных и начальных условиях. Зададим их в виде: T (Т, Xa ) = T1a (т), T- (т, Хь ) = T-ь (т), T (0, x) = Tw (x), T- (0, x) = T-o (x) . (6) Третье и четвертое уравнения неразрывности в (3) легко интегрируются с точностью до констант C1 (т), C2 (т) по длине потоков, откуда мы можем найти скорости хладагентов, зависящие от плотности потоков: ^т) = C ((Т) ))' = C (т) )) . (7) Р1 (T1 (Х,Т)) Р- (T2 (Х,Т)) Пространственные константы C1, C2 задаем через, как правило, известные плотности и скорости на входе или выходе потоков, например: C1 (т) = P1a (T1a (т)). (xa,т), C2 (т) = P-b (T2b (т))U2 (xb,т). Система (3) при условиях (5)-(7) принимает вид: 5-§=tV[kF(T- -T1)+q], ^-.2^=-1-[kF(T1 -T-)+q-]. ... [ w „ ^ A,-.IT=- kF(T1 -T-) + q- . (8) 5т 5x Cp1p^ 5т 5x Cp-р- Полученная математическая модель в виде уравнений (5)-(8), в отличие от описываемых ранее в литературе, позволяет для противоточных теплообменников моделировать температуру T(x^) с учетом не только переходных и фазовых превращений в хладагентах, но и пространственной распределенности физических параметров потоков, таких как плотность, скорость, теплоемкость и теплопередача. 3. Численное моделирование Для численного решения дифференциальных уравнений (8) для обеих сред использовалась одна и та же равномерная по пространству сетка Г xi = xa + / Ax, i = 0,1...n, Ax = x-a l с шагами по I n J Г TL времени . Для аппроксимации производных привлекалась неявная схе- V m J ма «прямоугольник» [11-13]. С й й дТ Согласно этой схеме значение производном по времени - усредняется по пространственным дт дТ точкам г и г + 1, а значение пространственной производном - усредняется по временным слоям j иj + 1. дх Кроме того, все коэффициенты уравнений определяются в центре рассматриваемого прямоугольника. Погрешность аппроксимаций такой схемы имеет второй порядок малости относительно Дт и Дх . Система (8) приобретает следующий конечно-разностный вид: (Tj+1 _tj + Tj+1 _tj )+(Tj+1 _Tj+1 + tj _tj) - 1 (tj+1/2 -tj+1/2)+q 2Дт\\ 1'i 1,г 1,i+1 11,i+1) +2Дх\\1>г+1 U U+1 11,г) ~ c p K ^ \\ 2,i+1/2 11,i+1/2) +41 S1 (9) P1 U А \\ 2,/ 2,/ 2,/+1 2,/ +1 1 r\\ А \\ 2,/+1 2,/ 2,/+1 2,/ / 2Дту ! 2Дх^ ! Cp2p2 Шаги аппроксимаций должны удовлетворять условию Куранта: Дх L_ (tj+1 _ tj + tj+1 _ tj ) V ( tj+1 _ tj+1 + tj _ tj ) - 1 kH (tj+1/2 _ tj+1/2 ) +q дД 2,/ 12,г+1 12,г+1) 2дД 2,/'+1 12,г + J2,г+1 12,i)~c p K ^ ^1,г+1/2 12,г+1/242 < Дт. там IV х,т (|и1( х,т)|,| и2 (х,т)|) Система уравнений (9) - это система 2(n + m) алгебраических уравнений, которую мы будем решать совместно для обеих сред в пространственных узлах сетки на каждом новом временном слое. Поскольку первая среда (в примере - воздух) движется слева направо, то его искомая температура должна находится в правой верхней точке конечно-разностного прямоугольника, т.е. в Т^-!. Соответственно для противоточной второй среды (в примере - кислород), движущейся справа налево, искомая температура должна находится в левой верхней точке Тj+1. Найти явным образом искомые температуры теплоносителей из уравнений (9) невозможно. Во-первых, потому что коэффициенты должны рассчитываться через неизвестную температуру Т+1/22 - (Т+1 + Т + Т--1 - Т+1)/4 . Во-вторых, как видно из правых частей уравнений (9), теплообмен определяется опять-таки неизвестной температурой противоточного теплоносителя. Для преодоления указанных проблем поступим следующим образом. Сначала выберем, какая из сред будет рассчитываться первой на каждом временном слое. Пусть это будет кислород, поскольку работа теплообменника заключается в охлаждении воздуха. Тогда коэффициенты второго уравнения (кислород) будем вычислять по трехточечному шаблону через (T2j++11,T2ji+1,T2/i), а для воздуха - через (T1J/+1,T1Ji,Т1]г+1) . Разности температур теплоносителей в правой части второго уравнения будем вычислять по данным только старого временного слоя, т.е. через (-Tj+1/2 _ T2ji+1/2 ) . Для первого уравнения (воздух), такой проблемы не возникает, так как температура кислорода будет уже известна. Теперь с учетом сделанных уточнений аппроксимаций можно из уравнений (9) выразить искомые температуры: j - Ц, -Ц^ _ + + М-, г - n _ 1...0, 1+ Y1 1+ Y1 1+ Y1 1+ Y1 1+ Y1 j -тг+11+Y1_0.5Y2 _T1r1_Y1+0.5Y2 + Tj 1+Y1+0.5Y2 + (10) 1+ Y1 + Y 2 ' 1+ Y1 + Y 2 ' 1+ Y1 + Y 2 + ^gl + 2Xq1 , / - 0...n _ 1, 1+ Y1 + Y 2 1+ Y1 + Y 2 иДт ДткИ Дт „ , ч „ , ч где y1 --, Y2 --, X --, S - первая среда (воздух), S2 - вторая среда (кислород). Дх S1,2Cp Р СрР 5. Результаты расчетов При расчете нестационарных процессов почти всегда возникает проблема определения начальных условий. Пусть в начальный момент времени противоточные среды не взаимодействуют т.е. k = 0 , тогда можно задать постоянные температуры, например T1 (0, x) = 315 К и T2 (0, x) = 90 К . На следующем временном слое j = 1, когда к ф 0 , начнется нестационарное течение. Рассмотрением именно такого нестационарного процесса мы в данной статье и ограничимся. В начальный момент времени примем давления в потоках p1 = 5-106Па, p2 = 2,4-106Па и константы С1 = 1285, С2 = 1717. Константа для воздуха задавалась при скорости и1 (xa ,т ) = 20 м/с, а для жидкого кислорода - при скорости и2 (xb ,т) = 1.5 м/с. Длина теплообменной поверхности аппарата принималась равной 9 м, площади поперечных сечений потоков S1 = 0,53м2, S2 = 0,18м2, площадь взаимодействия потоков F = 1100м2. Приведенные данные соответствуют типичному кожухотрубному теплообменному аппарату. Расчетная сетка задавалась Ax = 0,045 м и Ат = 0,00075 c, Т = 7 c - время полного установления процесса. Для указанного теплообменника воздух не меняет своего фазового состояния, ^ = 0, а вот кислород меняет. В этих условиях коэффициент теплопередачи между кислородом и воздухом нельзя считать постоянным по длине потоков. Будем задавать k = 6 -103Вт/м2К - для зоны жидкого кислорода, k = 4,5-103Вт/м2К - для зоны парожидкостной смеси кислорода, k = 3•103Вт/м2К - для зоны газообразного кислорода. На рис. 2 приведены результаты расчетов температуры воздуха (сплошная линия) и кислорода (пунктирная линия) в момент установления изначально нестационарного потока в рассматриваемом теплообменном аппарате. Две точки на температурных кривых для кислорода указывают на начало и конец парожидкостной зоны. В этой зоне, естественно, температура не меняется и равна Г2кип = 135 К. Рисунок 2, b получен при постоянном по всей длине потоков коэффициенте теплопередачи со средним значением k = 4,5 • 103 Вт/м2К. 350 Рис. 2, b. Температуры потоков при постоянном k Fig. 2, b. Flows temperatures at permanent k 50 300 . 250 200 150 100 \\ N \\ .........у N \\--- --- 0 50 100 150 Пространственная сетка, i 200 Рис. 2, а. Температуры потоков при переменном k Fig. 2, a. Flows temperatures at variable k Можно заметить, что распределенность коэффициента теплопередачи по длине теплообменно-го аппарата заметно влияет на процесс теплообмена. На рис. 3 показаны конечные распределения скоростей потоков. Видно, как кислород (пунктирная линия) ускоряется справа налево по мере кипения за счет изменения своей плотности от жидкости к газу. В то же время поток воздуха замедляется слева направо из-за существенного охлаждения и повышения своей плотности. Рис. 3, а. Скорости потоков при переменном k Fig. 3, a. Flows rate at a variable k Рис. 3, b. Скорости потоков при постоянном k Fig. 3, b. Flows rate at a permanent k Приведенные расчеты показывают, что учет фазовых переходов в хладагентах и учет изменения коэффициента теплопередачи вдоль теплообменного аппарата существенно влияют на физическую картину тепловых и гидродинамических процессов теплообменного аппарата. Заключение Полученная математическая модель (5)-(8) и ее конечно-разностное представление (10) позволяют осуществлять физически обоснованное численное моделирование одномерных нестационарных тепловых потоков в противоточных теплообменных аппаратах. Такая модель дает возможность учитывать распределенность физических параметров по длине потоков, таких как плотность, скорость, теплоемкость и теплопередача между потоками, и позволяет моделировать потоки противоточных хладагентов с фазовыми переходами.

Скачать электронную версию публикации