Лаб 9 МСУ. Контрольные вопросы что позволяют определить статические характеристики объекта Как осуществить линеаризацию статической характеристики Какие виды динамических характеристик вызнаете Какая взаимосвязь между импульсной и переходной характеристиками Какие методы аппроксимации динамических характеристи
44 КОНТРОЛЬНЫЕ ВОПРОСЫ Что позволяют определить статические характеристики объекта Как осуществить линеаризацию статической характеристики Какие виды динамических характеристик вызнаете Какая взаимосвязь между импульсной и переходной характеристиками Какие методы аппроксимации динамических характеристик вызнаете Чем отличается файл-функция от файла в пакете MATLAB? Какая функция может быть использована для проведения аппроксимации полиномом й степени Что осуществляет функция polyval в пакете MATLAB? Какая функция в пакете MATLAB может быть использована для поиска минимума функции нескольких переменных Лабораторная работа ПОСТРОЕНИЕ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ АНАЛИТИЧЕСКИМ МЕТОДОМ МОДЕЛИРОВАНИЕ ТЕПЛООБМЕННОЙ АППАРАТУРЫ Цель работы ознакомление с методами моделирования процесса теплообмена через стенку и с расчётом теплообменных аппаратов на ЭВМ. Задания работы Ознакомиться с принципами составления тепловых балансов потоков и разделяющей стенки в теплообменных и реакционных аппаратах химической технологии при теплообмене через стенку. Отметить особенности математических описаний при разных способах организации теплообмена. Выполнить вариант задания на ЭВМ. Содержание отчёта: 1. Вариант задания. 2. Математическая модель теплообменного аппарата. 3. Все расчётные формулы, используемые в работе. 4. Графики рассчитанных зависимостей для теплоносителей. 5. Тексты программ. КРАТКИЕ ТЕОРЕТИЧЕСКИЕ ПОЛОЖЕНИЯ Процесс передачи тепла через стенку весьма распространён в химической технологии и значительно влияет на протекание химических реакции во всех типах реакторов. Процесс передачи тепла в теплообменной 45 аппаратуре является основными служит для сообщения технологическому потоку нужной температуры. Выбирая различные способы оформления реакторов, можно влиять на интенсивность теплообмена между основным (реакционным) потоком и потоком хладагента или окружающей средой. При полном отсутствии теплообмена через стенку получают адиабатический реактор. Реакторы, имеющие теплообмен с внешней средой, относятся к политропическим. При рассмотрении процесса передачи тепла от одного теплоносителя к другому через стенку можно выделить несколько элементарных этапов переход тепла от горячего теплоносителя к более холодной стенке, поглощение тепла материалом стенки и её нагрев, распределение тепла по объёму стенки, переход тепла от стенки к холодному теплоносителю. Если процесс теплообмена протекает стационарно, то температура в каждой точке материала (теплоносителей и стенки) не изменяется во времени. Применение модели с сосредоточенными параметрами (те. когда пространственные координаты не входят в математическое описание) приводит к алгебраическим соотношениям между температурами в системе. Если, наоборот, температуры меняются во времени, математическое описание получается в виде системы обыкновенных дифференциальных уравнений (аргументом является время. Зависимость температур от геометрических координат обуславливает математическое описание статики в виде обыкновенных дифференциальных уравнений (если пространственная координата одна) или дифференциальных уравнений в частных производных. Независимыми переменными при этом являются пространственные координаты. Динамическая модель при наличии пространственно-распределённых эффектов описывается уравнениями в частных производных, причём одной из независимых переменных является время. Интенсивность перехода тепла от одного теплоносителя (например, горячего потока жидкости или газа) к другому (стенке) зависит от разности температур между ними, а также от теплового сопротивления. В рас- чётные уравнения, однако, обычно включают не сопротивление, а обратную величину α – коэффициент теплоотдачи, q – тепловой поток (ккал/ч или Вт) через поверхность площадью 1 м при разности температур (температурном напоре) 1 С. Полный тепловой поток q определяется произведением коэффициента теплоотдачи α , площади поверхности F и температурного напора ∆T: T F q ∆ α = . (4) Уравнение (4) применимо как к нагреванию стенки от горячей жидкости, таки, наоборот, к нагреванию холодной жидкости горячей стенкой при этом ∆T будет иметь разные знаки. 46 Если пренебречь распространением тепла в стенке, то теплопередачу от горячего потока жидкости к холодному, находящемуся по другую сторону стенки, можно представить как процесс преодоления тепловым потоком двух последовательных сопротивлений теплоотдачи – от горячего потока к стенке и от нагревшейся стенки к холодному потоку. Используя вместо сопротивлений коэффициенты теплоотдачи (α 1 и α 2 ), получаем выражение, определяющее коэффициент теплопередачи (K): 2 1 1 1 1 α + α = K (5) В практических расчётах часто используют коэффициент теплопередачи как характеристику интенсивности теплообмена между потоками T KF q ∆ = . (6) В тех случаях, когда коэффициенты теплоотдачи учитываются порознь, принимают усреднённую температуру стенки, разделяющей потоки. Иными словами, считают, что теплопроводность материала стенки настолько велика, что перепад температуры отсутствует. Коэффициенты теплоотдачи зависят от многих параметров, но наиболее сильно – от скорости потока, характера набегания жидкости на стенку, плотности и теплопроводности жидкости. При выполнении точных расчётов зависимость коэффициента теплоотдачи от параметров потока следует учитывать. Однако для большинства инженерных расчётов теплообменной аппаратуры и реакторов достаточны упрощённые представления. Для вывода уравнений математического описания процесса теплообмена через стенку следует рассмотреть тепловой баланс каждой среды, имеющей запас тепла. Он складывается из прихода и расхода тепла, которые определяют накопление тепла в объёме; накопление является временным процессом накопление = приход – расход. В статике ввиду равенства прихода и расхода тепла накопление равно нулю. Накопление связано с изменением температуры или для элементарного объёма SdTdl c p ρ , где ρ – плотность с – удельная теплоёмкость; V – объём; S – сечение потока d1 – элементарный участок потока. Приходи расход тепла могут определяться теплоотдачей (теплопередачей, а в случае проточной системы с распределёнными параметрами – притоком и уносом тепла с конвективным потоком. 47 Количество тепла, поступающее в аппарат с конвективным потоком, определяется как или для элементарного объёма за элементарное время τ d τ ρ Td c G p , где G – объёмный расход потока. Количество тепла, уходящее из рассматриваемого объёма с конвективным потоком, определяется следующим выражением τ ∆ + ρ ) ( T T c G p , или для элементарного объёма Приход тепла, определяемый теплопередачей τ − V T T V F K ) ( вн , или для элементарного объёма за элементарное время τ − dld T T V FS K ) ( вн , где Т вн – температура внешнего теплоносителя. С учётом полученных соотношений накопление тепла в системе составит τ − + τ ∆ + ρ − τ ρ = ∆ ρ ) ( ) ( вн T T KF T T c G T c G T V c p p p , или в элементарном объёме за элементарное время τ − + τ + ρ − τ ρ = ρ dld T T V FS K d dT T c G Td c G SdTdl c p p p ) ( ) ( вн Проведя несложные преобразования, получим уравнение теплового баланса, описывающее динамику теплообменников, во всём объёме которых происходит полное (идеальное) смешение частиц потока ) ( ) ( вн 0 T T KF T T c G d dT V c p p − + − ρ = τ ρ , (7) где Т, Т – температура потока на входе ив зоне идеального смешения. Соответственно для трубчатых теплообменников, работающих по принципу вытеснения, уравнение динамики будет выглядеть следующим образом V T T KF l T S c G T c p p ) ( вн − + ∂ ∂ ρ − = τ ∂ ∂ ρ (8) Ввиду того, что в статическом режиме накопление тепла в системе равно нулю, модель статики теплообменников смешения будет иметь вид 48 0 ) ( ) ( вн 0 = − + − ρ T T KF T T c G p , (9) статика трубчатых теплообменников описывается уравнением V c u T T KF l T p ρ − = ∂ ∂ ) ( вн , (10) где u – скорость движения потока. Пример 1. Теплообменник представляет собой тонкостенный змеевик, по которому в режиме идеального вытеснения движется охлаждаемый поток жидкости. Змеевик погружён вводу, непрерывно протекающую через сосуд, так что температура охлаждающей воды Т вн практически постоянна и равна 10 С во всём объёме. Требуется определить температуру на выходе потока, идущего по змеевику со скоростью u = 4 мс, если температура его на входе равна 95 С, длина трубки змеевикам, его сечением, коэффициент теплопередачи K = 1,16·10 4 Вт/(м 2 С. Теплоёмкость охлаждаемой жидкости c p = 2,93·10 3 Дж/(кг Се плотность ρ = 900 кг/м 3 . Параметры считать независящими от температуры изменение объёма не учитывать. Режим работы считать стационарным. Температура охлаждаемого потока Т подчиняется дифференциальному уравнению (10): S c u T T r K l T p ρ − π = ∂ ∂ ) ( 2 вн , (11) где l – длина r – радиус змеевика uS – объёмный расход потока. Начальное условие для уравнения (11) Т) = 95 С. Вычислим коэффициент уравнения 39 , 0 10 10 93 , 2 900 4 14 , 3 10 16 , 1 2 2 4 3 Тогда ) ( 39 , 0 вн T T l T − = ∂ ∂ (12) Уравнение (12) подлежит решению в пределах изменения независимой переменной 1 от 0 дом. Решение представлено на рис. 6. Температура на выходе потока составила С. Для численного интегрирования уравнения (12) в MATLAB необходимо создать два файла. В файле func1_T.m описывается функция для вычисления правой части дифференциального уравнения (12). В файле Tepl1.m осуществляются задание исходных данных для решения задачи и начального условия, а также решение дифференциального уравнения с использованием функции Matlab ode45. 49 0 0.4 0.8 1.2 1.6 2 45 55 65 75 85 95 T, o C l, м Рис. 6. Профиль температуры по длине теплообменника Функция ode45 в данном случае имеет следующий синтаксис [l, T] = ode45(@func1_T, [0 L], T0), где @func1_T – дескриптор функции правой части дифференциального уравнения [l, T] – матрица решений T, где каждая строка соответствует длине, возвращённой в векторе-столбце l; [0 L] – вектор, определяющий интервал интегрирования T0 – вектор начальных условий. Файл func1_T.m function dT = func1_T(l,T) % Функция правой части дифференциального уравнения Tvn=10; % Температура внешнего теплоносителя dT = 0.39*(Tvn-T); % Уравнение правой части % дифференциального уравнения Файл Tepl1.m L=2; % Длина трубки змеевика T0=95; % Начальное условие диф. уравнения [l,T] = ode45(@func1_T,[0 L],T0); % Функция решения дифференциального уравнения plot(l,T); % Построение зависимости T(l) Пример 2. Жидкость охлаждается в теплообменнике типа труба в трубе. Охлаждаемая жидкость и хладагент движутся параллельно (прямотоком. Требуется определить температуры потоков на выходе теплообменника, если начальная температура охлаждаемой жидкости равна 170 С, а хладагента 15 С. Плотность охлаждаемой жидкости и хладагента ρ = 900 кг/м 3 . Диаметры труб теплообменника внутренней D 1 = 0,1 м, наружной (для хладагента) D 2 = 0,3 м. Длина теплообменникам. Тепло- ёмкость жидкости и хладагента с = 3,35·10 3 Дж/(кг С. Объёмный расход охлаждаемой жидкости G 1 = 2,28·10 –4 мс, хладагента G 2 = 5,75·10 -4 мс, 50 коэффициент теплопередачи K = 4900 Вт/(м 2 ⋅ °С). Температурный профиль по длине для каждого из потоков определяется решением системы дифференциальных уравнений 1 1 1 1 2 1 1 ) ( p c G T T D K l T ρ − π = ∂ ∂ , (13) 2 2 2 2 1 1 2 ) ( p c G T T D K l T ρ − π = ∂ ∂ , (14) где T 1 и Т – температуры охлаждаемой и охлаждающей жидкости. Начальные условия Т) = 170 СТ С. После подстановки в уравнения (13) и (14) численных значений параметров получаем следующую систему ) ( 239 , 2 1 2 1 T T l T − = ∂ ∂ , (15) ) ( 888 , 0 2 1 2 T T l T − = ∂ ∂ (16) Графики решения системы уравнений математического описания статики теплообменника представлены на рис. 7. На нём изображены температурные профили вдоль теплообменника для обоих теплоносителей. Можно видеть, что движущая сила процесса сильно меняется по длине, поэтому эффективность использования различных участков теплообменника неодинакова. Температуры теплоносителей на выходе теплообменника равны T 1 (L) = 62,9 СТ С. Численное интегрирование уравнений (15) и (16) аналогично примеру. Создаются два файла func2_T.m – функция для вычисления правых частей дифференциальных уравнений, Tepl1.m – задание исходных данных, начальных условий и решение дифференциальных уравнений с использованием функции MATLAB ode45. 0 0.2 0.4 0.6 0.8 1 0 20 60 100 140 180 T 1 T 2 T, o C l, м Рис. 7. Изменение температур теплоносителей по длине прямоточного теплообменника 51 Файл func2_T.m function dT = func2_T(l,T) % Функция правых частей % дифференциальных уравнений dT = zeros(2,1); % Создание вектора температур dT(1)=2.239*(T(2)-T(1)); % Уравнение правой части го % диф. уравнения dT(2)=0.888*(T(1)-T(2)); % Уравнение правой части го % диф. уравнения Файл Tepl2.m L=1; % Длина теплообменника T1_0=170; % Начальное условие для го диф. уравнения T2_0=15; % Начальное условие для го диф. уравнения [l,T] = ode45(@func2_T,[0 L],[T1_0 T2_0]); Функция решения % дифференциальных уравнений plot(l,T(:,1),'r',l,T(:,2),'k'); % Построение зависимостей % T1(l) и Т) Пример 3. Смоделировать статический режим теплообменника типа труба в трубе, используя данные, приведённые в примере 2, для случая противотока. Принять полную длину теплообменникам. Тепловые процессы в противоточном теплообменнике подчиняются тем же закономерностям, что ив прямоточном. Поэтому математическое описание теплообменника записывается аналогично, однако формально однотипные уравнения для обоих теплоносителей имеют аргументы различного знака 1 1 1 1 2 1 1 ) ( p c G T T D K l T ρ − π = ∂ ∂ , (17) 2 2 2 2 1 1 2 ) ( ) ( p c G T T D K l T ρ − π = − ∂ ∂ (18) Существенным различием, отражающим иную организацию потоков теплоносителей, является принципиально другое задание условий решения уравнений (17) и (18) по сравнению с заданием при решении уравнений) и (14). Совместное интегрирование уравнений (17) и (18) возможно лишь водном направлении либо применяющемся от 0 до L, либо в обратном – от L до 0. При этом в любом случае оговорено лишь одно начальное условие, второе остаётся неизвестным. Известно лишь, к какому значению в конце решения должна подойти вторая переменная. Для решения задачи воспользуемся последним обстоятельством и попытаемся отыскать неизвестное начальное условие Т) с таким расчё- том, чтобы условие, заданное для конца решения (граничное условие, было выполнено, те. T 2 (L) = 15 С. Такие задачи при малом числе условий, подлежащих определению, обычно решают методом проб и ошибок. 52 Задачей поиска начального условия Т) является выполнение граничного условия T 2 (L) при интегрировании системы уравнений (17) и (18). Рисунок 8, а иллюстрирует процесс поиска неизвестного начального условия Т. Кривая 1 отражает профиль температуры Т, полученный в предположении, что хладагент нагреется до 57 С (значение взято произвольно. Видно, что это предположение неудачно, так как не получено ожидаемое значение T 2 (L) = 15 С. 0 0 .5 1 1 .5 2 2 .5 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 1 0 0 T 2 , o C l, мамб) Рис. 8. Результаты моделирования на ЭВМ противоточного теплообменника а – серия кривых решения системы уравнений (17) и (18), полученная в процессе определения недостающего начального условия б – изменение температур теплоносителей по длине теплообменника 53 Приняв Т) = 65 С, получаем кривую 2. Она также неудовлетворительна, так как вместо Т) = 15 С получаем лишь 4 С. Кривые 3 и 4 также не соответствуют искомому начальному условию в первом случае входная температура T 2 (L) получилась заниженной, во втором – завышенной. По-видимому, искомое начальное условие лежит между двумя последними проверенными значениями (70 и 80 С. Тщательно исследуя намеченный диапазон задания начального условия, находим Т) = 75 С. На рисунке 8, б показан результат решения задачи. Выходная температура охлаждаемого потока T 1 (L) достигает 18 С. Численное интегрирование уравнений (17) и (18) аналогично примеру. Файл func3_T.m function dT = func3_T(l,T) % Функция правых частей % дифференциальных уравнений global b1 b2; % Описание глобальных переменных dT = zeros(2,1); % Создание вектора выходных координат dT(1)=b1*(T(2)-T(1)); % Уравнение правой части го % диф. уравнения dT(2)=-b2*(T(1)-T(2)); % Уравнение правой части го % диф. уравнения Файл Tepl3.m global b1 b2; % Описание глобальных переменных L=2.5; % Длина теплообменника T1_0=170; % Начальное условие для го диф. уравнения T2_0=80; % Начальное условие для го диф. Уравнения ro=900; % Плотность охлаждаемой жидкости и хладагента D1=0.1; % Диаметры внутренней трубы теплообменника D2=0.3; % Диаметры наружной трубы теплообменника cp=3.35e3; % Теплоемкость жидкости и хладагента G1=2.28e-4; % Объемный расход охлаждаемой жидкости G2=5.75e-4; % Объемный расход хладагента K=4900; % Коэффициент теплопередачи % Расчет коэффициентов модели b1=K*pi*D1/(ro*cp*G1); b2=K*pi*D1/(ro*cp*G2); [l,T] = ode45(@func3_T,[0 L],[T1_0 T2_0]); Функция решения % дифференциальных уравнений plot(l,T(:,1),'b',l,T(:,2),'k'); % Построение зависимостей % T1(l) и Т) Пример 4. Смоделировать переходный режим теплообменника типа смешение – смешение. Теплообменник представляет собой двухкамерную ёмкость. В первую камеру ёмкости поступает охлаждаемая жидкость, во вторую – хладагент. Для обеспечения однородного распределения температуры по объё- му в камерах установлены мешалки. Плотность охлаждаемой жидкости 850 кг/м 3 , а хладагента 920 кг/м 3 . Объёмы камер равны и составляют 2,5 м каждая. Объёмный расход теплоносителя 4,12·10 –3 мс, хладагента 5,43·10 –3 мс. Теплоёмкости жидкости и хладагента соответственно 3,75·10 3 Дж /(кг ⋅ °С) и 3,14·10 3 Дж/(кг С. Поверхность теплообмена составляет 4 м, а коэффициент теплопередачи равен K = 4360 Вт/(м 2 С. Температура охлаждаемой жидкости на входе меняется скачкообразно от 115 до 200 С, а температура хладагента от 10 до 15 С. Изменение температуры повремени для каждого потока определяется решением системы ) ( ) ( 1 2 1 0 1 1 1 1 1 1 1 1 T T KF T T c G d dT V c p p − + − ρ = τ ρ , (19) ) ( ) ( 2 1 2 0 2 2 2 2 2 2 2 2 T T KF T T c G d dT V c p p − + − ρ = τ ρ , (20) где 0 1 T и 0 2 T – температуры жидкости и хладагента на входе в теплообменник. Система уравнений (19) и (20) решается при начальных условиях * 1 1 ) 0 ( T T = ; * 2 Для определения * 1 T и * 2 T , соответствующих номинальному статическому режиму, составляются уравнения теплового баланса стационарного режима работы теплообменника 0 ) ( ) ( 1 2 1 0 1 1 1 1 = − + − ρ T T KF T T c G p , 0 ) ( ) ( 2 1 2 0 2 2 2 После подстановки численных значений получаем следующую систему, которая может быть решена численно или аналитически относительно любых двух переменных, входящих в эти уравнения. Температуры теплоносителей на выходе теплообменника в установившемся состоянии составили * 1 T = 74,46 С, * 2 T = 43,93 С. Таким образом, получим следующую систему 55 ) ( ) ( 1 2 1 1 1 1 0 1 1 1 1 T T V c KF T T V G d dT p − ρ + − = τ , (21) ) ( ) ( 2 1 2 2 2 2 0 2 1 2 2 T T V c KF T T V G d dT p − ρ + − = τ (22) Начальные условия Т) = 74,46 СТ С. Графики решения системы уравнений математического описания динамики теплообменника представлены на рис. 9. На нём изображены изменения температур во времени для обоих теплоносителей на выходе теплообменника. Температуры теплоносителей на выходе теплообменника установились через 2000 си составили T 1 = 126,6 СТ С. Численное интегрирование уравнений (21) и (22) аналогично предыдущим примерам 1 – 3. Отличие состоит в том, что интегрирование уравнений с помощью функции ode45 производится повремени, с T 1 T 2 Рис. 9. Изменение температур теплоносителей во времени теплообменника типа смешение – смешение Файл func4_T.m function dT = func4_T(time,T) % Функция правых частей % диф. уравнений global a1 b1 a2 b2 T1_vx T2_vx; % Описание глобальных переменных dT = zeros(2,1); % Создание вектора температур dT(1)=a1*(T1_vx-T(1))+b1*(T(2)-T(1)); % Уравнение правой % части го диф. уравнения 56 dT(2)=a2*(T2_vx-T(2))+b2*(T(1)-T(2)); % Уравнение правой % части го диф. уравнения Файл Tepl4.m global a1 b1 a2 b2 T1_vx T2_vx; % Описание глобальных % переменных tk=4000; % Время интегрирования T1_0=74.46; % Начальное условие для го диф. уравнения T2_0=43.93; % Начальное условие для го диф. уравнения T1_vx=200; % Температура охлаждаемой жидкости на входе T2_vx=10; % Температура хладагента на входе ro1=850; % Плотность охлаждаемой жидкости ro2=920; % Плотность охлаждаемой хладагента cp1=3.75e3; % Теплоёмкость охлаждаемой жидкости cp2=3.14e3; % Теплоёмкость хладагента G1=4.12e-3; % Объёмный расход охлаждаемой жидкости G2=5.43e-3; % Объёмный расход хладагента K=4360; % Коэффициент теплопередачи F=4; % Площадь поверхности теплообмена V=2.5; % Объём камер % Расчёт коэффициентов модели a1=G1/V; a2=G2/V; b1=K*F/(ro1*cp1*V); b2=K*F/(ro2*cp2*V); [time,T] = ode45(@func4_T,[0 tk],[T1_0 T2_0]); Функция % решения дифференциальных уравнений plot(time,T(:,1),'b',time,T(:,2),'k'); % Построение % зависимостей T1(time) и Т) ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ 1. Составить математическое описание теплообменного аппарата варианты теплообменных аппаратов и исходные данные приведены в табл. 6, 7, 8). 2. В зависимости от варианта смоделировать на ЭВМ статический и или) динамический режимы теплообменника и определить температурные зависимости для всех теплоносителей. КОНТРОЛЬНЫЕ ВОПРОСЫ Чем определяется интенсивность перехода тепла от одного теплоносителя к другому |