А. И. Данилов Московский государственный университет пищевых производств Предисловие Данный практикум
| Вид материала | Практикум |
- «Московский государственный университет пищевых производств», 995.39kb.
- Совершенствование качества системы управления персоналом на предприятиях общественного, 565.85kb.
- Управление затратами предприятий пищевых производств в условиях формирования вертикально-интегрированных, 253.14kb.
- «Московский государственный университет пищевых производств», 382.12kb.
- Московский Государственный Университет пищевых производств Т. М. Панченко методические, 26.11kb.
- Комплексная технология переработки гречихи с утилизацией лузги, 323.59kb.
- Разработка технологических решений по совершенствованию сортовых помолов ржи, 327.16kb.
- Разработка системы обеспечения сенсорного качества соковой продукции, 392.87kb.
- Совершенствование технологии кристаллического сахара повышенной чистоты, сохраняемости, 360.85kb.
- Московский Государственный Университет пищевых производств Косикова Ю. А. методические, 186.42kb.
Работа 4. Моделирование объектов с распределёнными параметрами
Содержание
- Введение
- Математическая модель процесса
- Моделирование процесса при неизменной скорости
- Моделирование процесса при возмущении по скорости
- Задание и контрольные вопросы
- Дробно-рациональная аппроксимация трансцендентной передаточной функции
- Задание
1. Введение
В химической, нефтяной, газовой, пищевой и других отраслях промышленности имеется обширный класс объектов регулирования, процессы в которых характеризуются относительно большим временем пребывания потоков в технологическом аппарате и в гидродинамическом отношении близких к аппаратам идеального вытеснения.
Большая группа процессов, происходящих в указанных массообменных и теплообменных аппаратах, характеризуется постоянством внешнего воздействия по направлению движения обрабатываемого потока. Сюда относятся процессы сушки в ленточных сушилках с поперечным обдувом воздухом движущегося потока материала; процессы абсорбции при малых концентрациях абсорбируемого компонента или плохой растворимости его в абсорбенте; химические реакции первого порядка, протекающие в близких к изотермическим трубчатых реакторах; процессы теплообмена в аппаратах с постоянной по длине температурой греющего (охлаждающего) теплоносителя и т. п.
На этапе разработки систем автоматического регулирования таких объектов, часто требуется осуществлять на ЭВМ математическое моделирование их динамических режимов. Если на объект воздействуют сравнительно малые возмущения, то динамические режимы можно описывать линейными дифференциальными уравнениями в частных производных.
Известно несколько методов математического моделирования таких динамических режимов. Излагаемый далее способ моделирования основан на анализе структуры трансцендентных передаточных функций. Этот способ позволяет исследовать динамику объекта с распределенными координатами по всем каналам, используя практически одни и те же динамические звенья; при этом важно то, что получение динамической модели для различных длин аппарата сводится к нетрудоемкой операции изменения коэффициентов усиления одних и тех же динамических звеньев модели.
^ 2. Математическая модель процесса
Так как законы тепло- и массопередачи принципиально идентичны, то в дальнейшем все рассуждения будут иллюстрированы на примере процесса теплообмена.
Рассмотрим теплообменники (рис. 1):
а) с постоянной по длине температурой греющего агента
θ3(x) = const;
б) с постоянным тепловыделением по длине теплообменника
q(x) = const.
Примем следующие допущения:
- Свойства обрабатываемого потока неизменны во времени и по длине аппарата.
- Температура стенки берется как среднеинтегральная по толщине.
- Тепловой поток в аксиальном направлении пренебрежимо мал.
- Коэффициенты теплоотдачи постоянны по длине аппарата.
Выделим на расстоянии х от входа обрабатываемого потока в аппарат элемент dx (рис.1а). Составим дифференциальный тепловой баланс для элементарного кольца стенки за время dt:
Аналогично составляем тепловой баланс для элемента обрабатываемого потока внутри трубы:
Разделим левую и правую части уравнений на dt. Так как поток движется, то полная производная
будет равна
Проведя сокращение и учитывая выражение для производной
, получим систему двух дифференциальных уравнений в частных производных (1), (2) с краевыми и начальными условиями (3). Необходимо самостоятельно проделать все преобразования.  |  | (1) |
 | (2) | |
 | (3) |
Аналогично можно составить тепловой баланс в дифференциальной форме для теплообменника с внутренним источником тепла (рис. 1б):
| |  | (4) |
 | (5) | |
 | (6) |
Здесь f 1(t), f 3(t) и
- произвольные функции времени.Преобразуем системы (1)-(3) и (4)-(6) по Лапласу относительно t с учетом начальных условий
| |  | (7) |
 | (8) | |
 | (9) |
| |  | (10) |
 | (11) | |
| | (12) |
Так как системы (7)-(9) и (10)-(12) идентичны, то далее будем оперировать только с первой из них.
Найдем решение системы (7)-(9) относительно температуры θ1 обрабатываемого потока, принимая во внимание, что температура греющего агента θ3=f3(t) есть величина, не зависящая от координаты x и изменяющаяся только во времени. С краевым условием (9), получаем решение системы (7)-(9)
 | (13) |
 | |
 |
или в более короткой записи
 | (14) |
или
 | (15) |
Необходимо самостоятельно уметь проделывать преобразования, переводящие систему типа (1)-(3) в систему (7)-(9) . Необходимо также уметь самостоятельно получать решение (13) дифференциального уравнения типа (7)-(9), где независимая переменная x - текущая длина аппарата.
Рис. 1
^ 3. Моделирование процесса при неизменной скорости
На основании (15) можно составить блок-схему моделирования на Simulink΄e подобного рода процессов, рис. 2.
Реализация звена с передаточной функцией W1=1/M(p) не вызывает затруднений. Рассмотрим возможность реализации звена
Учитывая, что порядок многочлена ^ М(р) выше порядка полинома N(p), и проведя деление M(p)/N(p), получим
 | (16) |
Таким образом, можно записать
 | (17) |
где
- чистое запаздывание, а
- постоянный коэффициент. Остается реализовать звено с передаточной функцией
 | (18) |
или
, которая входит сомножителем в (17). Это нетрудно выполнить, разложив правую часть выражения (18) в ряд  | (19) |
где
Рис. 2
На примере, рис. 3, представлена структурная схема, позволяющая исследовать с помощью Simulink-модели изменение температуры θ1(t) в зависимости от возмущений f1(t), f3(t).
При моделировании были взяты следующие значения постоянных времени, T12=T21=T23=1 и время запаздывания
с. Выражение (15) при этом примет вид
 | (20) |
Для реализации звена с трансцендентной передаточной функцией
 | (21) |
понадобилось три интегратора. Статическая ошибка, если ограничиться первыми четырьмя членами разложения для (21), составит
что вполне допустимо.
Ранее было указано, что количество членов разложения зависит от величины безразмерного комплекса β. Это ясно из анализа разложения (19) при p->0, т.е. при t =
, когда наступает установившееся (стационарное) состояние. Поскольку каждый член (19) представляет собой апериодическое звено n-го порядка (не колебательное!!!), смело можно утверждать, что погрешность "в динамике" от отбрасывания членов ряда, начиная с "n+1" члена, не будет превосходить погрешности "в статике", т.е.:
Указанная формула позволяет для заданной погрешности в статике
определить число необходимых членов разложения в (19), т.е. количество необходимых интеграторов, рис. 3. Для величин β>1 сходимость ряда ухудшается и может понадобиться очень большое (неприемлемое) количество интеграторов. Однако, по результатам проведенных научных исследований, можно высокоточно аппроксимировать звено
дробно-рациональными передаточными функциями, см. разд. 6.
^ 4. Моделирование процесса при возмущении по скорости
В большинстве случаев важно определить влияние на θ1 изменение скорости потока
. Эта величина входит в уравнение (1) в качестве параметра. Пусть скорость потока была 0 и соответственно температуры θ10 и θ20. Изменим скорость так, чтобы она стала 
, где δ(t) достаточно малая, изменяющаяся во времени, величина. Тогда, предполагая, что влияние возмущения на решение мало, можно записать уравнение (1), (2) в виде:  |  | (22) |
 | (23) |
Пренебрегая членами второго порядка малости, систему (22), (23) (с учетом начальных и краевых условий) можно записать в виде двух систем дифференциальных уравнений:
исходной (невозмущенной)
 |  | (24) |
 | (25) | |
 | (26) |
и системы в отклонениях
 |  | (27) |
 | (28) | |
 | (29) |
Производная по x решения системы (24-26) входит сомножителем в возмущающую функцию уравнения (27). Определим ∆θ1 из системы (27-29). При этом будем предполагать, что функции f1(t) и f3(t) есть постоянные во времени соответственно f10 и f30 и, что переходные процессы в аппарате закончились до нанесения возмущения δ
(t) по скорости. Тогда из уравнения (15) имеем
 | (30) |
Подставляя (30) в (27) и решая систему (27-29) таким же образом, как систему (1-3), получим операторное уравнение
 | (31) |
Моделирование этого операторного соотношения аналогично моделированию соотношения (15). Выражение
есть постоянное число; δ/0 - относительное возмущение по скорости. Но, внимание, в формулу (31) входит множитель 1/р и на первый взгляд может показаться, что после нанесения возмущения по скорости, величина ∆θ1(t) будет бесконечно возрастать.
Однако при тщательном рассмотрении видно, что
Раскрыв неопределенность, можно показать, что соотношение (31) имеет конечный предел при р->0
На рис. 4 представлена схема моделирования передаточной функции по изменению скорости потока. Особенностью этой схемы является необходимость тщательной установки коэффициента
так, чтобы он в точности был равен коэффициенту усиления звена 
в статике.Выражение (31) получено без учета изменения коэффициента теплоотдачи при изменении скорости. Однако можно показать, что выражение для ∆θ1(p,x), в случае учета изменения коэффициента теплоотдачи, будет идентичным
 | (32) |
Статический коэффициент усиления при этом будет меньше
Это можно было предположить, так как недогрев за счет увеличения скорости потока компенсируется увеличением коэффициента теплоотдачи.
Результаты моделирования рассмотренного ранее примера приведены на рис. 5, где показаны:
- реакции системы (20) на единичные скачкообразные возмущения по f3 и f1 (рис.5а);
- реакции системы (31) и (32) при ступенчатом изменении скорости потока (рис. 5б).
Постоянные времени имели значения:
Рис. 5. Кривые разгона: 1 - по каналу f1;
2 - по каналу f3; 3 - по скорости δ
для (31);4 - по скорости δ:
для (32)При выводе математических моделей были приняты следующие обозначения:
Координаты и параметры
x- текущая координата аппарата, м;
l - длина аппарата, м;
t - текущее время, с;
D,d - диаметр, м;
с - теплоемкость,
γ - плотность,
α - коэффициент теплоотдачи,
- скорость движения,0 - постоянная начальная скорость,δ
- возмущение скорости,θ - температура,
q - удельный тепловой поток,
k - коэффициент,
f, ψ - возмущающие функции;
- преобразование Лапласа по переменной t соответственно для θ, f , δ;δθ - отклонение температуры от номинального значения вследствие изменения скорости δ
;^ БПЗ - блок постоянного запаздывания.
Индексы
1 - обрабатываемый поток; 2 - стенка; 3 - греющий теплоноситель; 10, 20, 30 - номинальное значение параметра соответствующей среды; 12 - обрабатываемый поток - стенка;
21 - стенка - обрабатываемый поток; 23 - стенка - греющий теплоноситель.
Постоянные времени
Многочлены от р
^ 5. Задание и контрольные вопросы
Составить математическую модель объекта с распределенными параметрами рассмотренного типа (скорость неизменна). Постоянные времени T12, T21, T23, длина аппарата x и скорость движения обрабатываемого потока
заданы в приложении.Считать условно, что температура пара f3 может меняться в пределах от 0 до 100˚C, а входная температура f1 - от 0˚C до 10˚C
В пояснительной записке следует привести:
а) операторное соотношение (20) для соответствующего варианта;
б) блок-схему моделирования;
в) структурную схему модели;
г) результаты моделирования: кривые разгона θ1 по каналам f1 и f3.
Знать ответы на следующие контрольные вопросы:
- Как определить время запаздывания, которое необходимо реализовать, зная величины x и ?
- При каком значении безразмерной величины
коэффициент усиления передаточной функции W2 будет меньше 0,01, т. е. изменение входной температуры f1 практически не будет влиять на температуру выхода θ1 ?
- Каким образом будут изменяться коэффициенты усиления по каналам f1 и f3 при изменении скорости обрабатываемого потока (например, от 0.01 до 10 м/с)?
- Сколько членов разложения в ряд Тейлора передаточной функции W3 понадобится, чтобы ошибка в коэффициенте усиления не превысила % (в статике)?
- Как подсчитать коэффициенты модели a, b, m, c, d, e, f, g, h, рис.3, имея в качестве исходных данных постоянные времени T12, T21, T23, T.
- Как подсчитать коэффициенты модели a, b, m, c, d, e, f, g, h, рис.3, имея в качестве исходных данных c1, γ1, d, α21, c2, γ2, D, α32, x, .
^ 6. Дробно-рациональная аппроксимация трансцендентной передаточной функции
Введем безразмерную частоту
и оператор Р=α-р. Тогда 
будет передаточной функцией от одного параметра β. Сделаем дополнительно нормирующее преобразование  | (33) |
которое обеспечит изменение
от 1 до 0 при изменении частоты ω от 0 до .Далее приводятся данные, позволяющие с высокой точностью аппроксимировать трансцендентный дифференциальный оператор
дробно-рациональным 
для любого значения β.Ниже указаны четыре типа передаточной функции
, каждый из которых имеет наибольшую точность аппроксимации в своем числовом диапазоне параметра β.^ Внимание !!! Чтобы перейти к исходной передаточной функции W3 от ее высоко точной аппроксимации
необходимо проделать обратное нормирование, т.е.  | (34) |
| Для   где  |
| Для   где    |
| Для   где     |
| Для   где A1 = ^ 2; C = (целая часть числа β/2); A2 = 1 - (дробная часть числа β/2). |
Задание
Табл. заданий
| № вар | Т12 | Т21 | Т23 | | X |
| 1 | 8.4 с | 7.00 с | 1.33 с | 0.4 м/с | 4.8 м |
| 2 | 3.2 мин | 6.35 мин | 1.15 мин | 0.05 м/мин | 0.5 м |
| 3 | 7.1 мин | 5.10 мин | 0.63 мин | 1 м/мин | 8 м |
| 4 | 5.5 мин | 1.90 мин | 0.34 мин | 1.2 м/мин | 7.2 м |
| 5 | 4.6 с | 2.80 с | 0.24 с | 1.6 см/с | 6.4 см |
| 6 | 2.7 мин | 1.08 мин | 0.12 мин | 30 м/мин | 60 м |
| 7 | 0.56 ч | 0.48 ч | 0.21 ч | 200 см/ч | 210 см |
| 8 | 0.60 мин | 0.24 мин | 0.04 мин | 16 м/мин | 12.8 м |
| 9 | 0.23 с | 0.38 с | 0.62 с | 24 см/с | 14.4 см |
| 10 | 0.06 с | 0.178 с | 1.15 с | 1.0 мм/с | 0.4 мм |
| 11 | 0.11 мин | 0.244 мин | 0.27 мин | 30 м/мин | 6.0 м |
| 12 | 0.09 ч | 0.127 ч | 0.07 ч | 1115 см/ч | 100 см |