«МОДЕЛИРОВАНИЕ СИСТЕМ Издательство ФГБОУ ВПО ТГТУ Учебное издание ЕЛИЗАРОВ Игорь Александрович, МАРТЕМЬЯНОВ Юрий Федорович, СХИРТЛАДЗЕ Александр Георгиевич, ТРЕТЬЯКОВ Александр Александрович МОДЕЛИРОВАНИЕ СИСТЕМ Учебное ...»

Порядок реакции или стадии будем обозначать цифрой над (или под) стрелкой, указывающей направление реакции; здесь же будем приводить обозначение константы скорости. Так, запись обозначает реакцию, имеющую 1-й порядок по веществу А и В и константу скорости k.

На практике встречаются реакции самых разнообразных порядков:

целочисленных и дробных, нулевого, изредка и отрицательных.

Пример 5. Скорость реакции по веществам.

Рассмотрим реакцию имеющую первый порядок по реагенту А и константу скорости k. Для этой реакции, в соответствии с формулой (31), имеем Сложная реакция состоит из ряда стадий, каждая из них имеет свою скорость, и нельзя достаточно точно определить скорость всей сложной реакции в целом.

В то же время любое из участвующих в реакции веществ образуется или расходуется с определённой скоростью, причём знание этой скорости обязательно при расчёте процесса.

Скорость сложной реакции по любому веществу i равна алгебраической сумме скоростей всех стадий по этому веществу (с учётом стехиометрических коэффициентов):

Здесь j – номер стадии; m – число стадий.

Скорость стадии по веществу определяется формулой (31). Если вещество в стадии не участвует, его стехиометрический коэффициент в этой стадии равен нулю.

Пример 6. Скорости сложной реакции по веществам.

Запишем кинетические уравнения для веществ, участвующих в реакции:

Вещество А участвует в 1-й и 2-й стадиях со стехиометрическими коэффициентами –1 и –1. Поэтому Соответственно:

Если химическая реакция (простая или сложная) протекает в замкнутой системе при постоянном объёме, то ход её описывается системой дифференциальных уравнений. Действительно, совместное рассмотрение уравнений (23), (33) и (36) приводит к формуле:

причём число таких уравнений равно числу веществ, участвующих в реакции (реагентов и продуктов). При этом для сложных химических реакций число уравнений может стать значительным, что существенно затруднит возможность использования математической модели при исследовании такого процесса.

Число уравнений математической модели, описывающих изменение количеств компонентов химической реакции с течением процесса, можно сократить до некоторого минимума путём написания указанных уравнений только для ключевых компонентов реакции. Поведение остальных компонентов реакции, не вошедших в число ключевых, при этом описывается простыми стехиометрическими соотношениями, представленными через количества ключевых компонентов.

Для иллюстрации этих правил воспользуемся конкретным примером следующей химической реакции:

где,,, – стехиометрические коэффициенты.

Кинетика процесса полностью описывается, если заданы скорости r1, r 2, r 3 и r 4 всех стадий. В этом случае скорость образования любого компонента, участвующего в сложной реакции, можно записать через скорости стадий, в которых участвует данный компонент, и соответствующие стехиометрические коэффициенты с учётом правила знаков:

Составим матрицу стехиометрических коэффициентов для системы уравнений (5.17):

где коэффициент, стоящий на пересечении i-й строки и k-гo столбца, представляет собой стехиометрический коэффициент i-гo компонента в k-й стадии реакции, взятый с учётом знака.

Задача выбора ключевых компонентов сводится к нахождению ранга матрицы стехиометрических коэффициентов (40). Под рангом матрицы понимается порядок наименьшего определителя, который можно построить из матрицы. Ранг матрицы и, следовательно, число ключевых компонентов можно определить известными методами матричного исчисления.

Для случая матрицы (40) можно поступить следующим образом. Будем искать ранг матрицы, последовательно увеличивая порядок определителей. Очевидно, что определитель первого порядка, т.е. содержащий только один элемент, можно построить легко, для чего необходимо взять любой элемент матрицы, отличный от нуля.

Теперь переходим к отысканию определителя второго порядка, отличного от нуля. Берём первый определитель, стоящий в левом верхнем углу матрицы (40):

Он оказывается равным нулю. Нетрудно проверить, что и любой определитель, составленный из строк первых двух столбцов матрицы, также равен нулю, что в данном случае указывает на линейную зависимость между первым и вторым столбцами матрицы (40). Действительно, второй столбец может быть получен из первого простым умножением на –1.

Аналогично обстоит дело с определителями второго порядка, построенными из строк третьего и четвёртого столбцов.

Остаётся рассмотреть определители второго порядка, которые могут быть построены из строк второго и третьего столбцов. Легко убедиться, что для этих столбцов любая комбинация двух строк даёт определители, отличные от нуля, например определители:

Итак, получено, что ранг матрицы (40) равен 2 и, следовательно, в качестве ключевых компонента можно использовать два компонента.

Пусть в качестве ключевых компонентов выбраны А и D; тогда, как нетрудно проверить, скорости образования компонентов В и С могут быть выражены через скорости образования rA и rD :

Если же в качестве ключевых компонентов выбрать А и В, то для rC и rD соответственно получим:

Аналогично можно найти соотношения и для любой другой выбранной пары ключевых компонентов: А и С; В и С; В и D; С и D.

Уравнения, описывающие изменение концентраций компонентов со временем, имеют в данном случае вид:

С учётом соотношений (43) второе и третье уравнения системы (45) могут быть заменены выражениями:

после интегрирования которых получим:

где CA0, CB0, CC0 и CD0 – начальные концентрации компонентов.

Отсюда окончательно находим соотношения:

которые могут использоваться в системе уравнений математического описания химической реакции вместо второго и третьего дифференциального уравнения системы (45).

Пример 7. Моделирование химической реакции омыления этилсокцената едким натром.

Для упрощения записи введём обозначения для веществ, участвующих в реакции:

С2Н5ООС – СН2 – СН2 – СООС2Н5 – A;

C2H5OOC – CH2 – CH2 – COONa – С;

NaOOC – СН2 – СН2 – COONa – E.

Тогда схему реакции можно записать так:

Реакция протекает в жидкой фазе, т.е. практически изохорно. При 0 °С k1 = 2,24 м3/(кмольс); k2 = 0,33 м3/(кмольс).

Вследствие изохорности реакции можно применить уравнение (26):

По схеме скорости реакции равны Нетрудно убедиться, что ранг матрицы стехиометрических коэффициентов равен 2, вследствие чего выберем в качестве ключевых любые два компонента (например, компоненты А и Е). Тогда скорости остальных компонентов можно выразить следующими соотношениями:

С учётом последних соотношений математическое описание реакции будет выглядеть следующим образом:

На рисунке 10 приведены графические результаты моделирования – полученные кинетические кривые для начальных условий: СА0 = 1 кмоль/м3;

СB0 = 0,5 кмоль/м3.

Необходимо отметить, что без применения вычислительной техники расчёт очень громоздок; использование ЭВМ позволяет весьма просто и быстро рассчитать кинетические кривые для любых начальных концентраций и констант скоростей.

Рис. 10. Графические результаты моделирования реакции омыления Программа на Matlab, решающая систему уравнений (5.25) представлена ниже.

dt=0.001;

time=[0:dt:3];

n=length(time);

% Создание массивов концентраций Ca=zeros(1,n);

Cb=zeros(1,n);

Cc=zeros(1,n);

Cd=zeros(1,n);

Ce=zeros(1,n);

% Начальные условия Ca(1)=1;

Cb(1)=0.5;

Cc(1)=0;

Cd(1)=0;

Ce(1)=0;

% Константы скоростей реакции k1=2.24; k2=0.33;

% Интегрирование дифференциальных уравнений методом Эйлера for i=2:n Ca(i)=Ca(i-1)+dt*(-k1*Ca(i-1)*Cb(i-1));

Ce(i)=Ce(i-1)+dt*(k2*Cc(i-1)*Cb(i-1));

Cb(i)=Cb(1)+(Ca(i-1)-Ca(1))-(Ce(i-1)-Ce(1));

Cc(i)=Cc(1)-(Ca(i-1)-Ca(1))-(Ce(i-1)-Ce(1));

Cd(i)=Cd(1)-(Ca(i-1)-Ca(1))+(Ce(i-1)-Ce(1));

end % Построение кинетических кривых plot(time,Ca,'k', time,Cb,'b', time,Cc,'r', time,Cd,'g', time,Ce,'m');

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

1. Выразить уравнения скоростей реакции по каждому веществу, участвующему в процессе (варианты контрольных заданий и исходные данные приведены в табл. 9, 10).

2. Провести выбор ключевых компонентов реакции.

3. Составить математическое описание химической реакции с учётом пункта 2.

4. Определить кинетические зависимости для веществ, участвующих в процессе.

CA0, кмоль/м3 CB0, моль/м

КОНТРОЛЬНЫЕ ВОПРОСЫ

1. Что такое скорость химической реакции?

2. Кинетика простых реакций.

3. Кинетика сложных реакций.

4. Каким образом можно сократить число уравнений математической модели, описывающих изменение количеств компонентов химической реакции с течением процесса?

5. Как определяется порядок химической реакции?

6. Объясните физический смысл уравнения Аррениуса.

МОДЕЛИРОВАНИЕ ПРОСТЫХ ГИДРАВЛИЧЕСКИХ СИСТЕМ

Цель работы: освоить методику математического моделирования гидравлических систем на ЭВМ.

Задания работы:

1. Ознакомиться с основными принципами моделирования гидравлических систем.

2. Выполнить вариант задания.

Содержание отчёта:

1. Вариант задания.

2. Математическая модель гидравлической системы.

3. Все расчётные формулы, используемые в работе.

4. Тексты программ и результаты вычислений.

КРАТКИЕ ТЕОРЕТИЧЕСКИЕ ПОЛОЖЕНИЯ

К простым гидравлическим системам (рис. 11) относятся технологические схемы трубопроводов, для которых принимаются следующие допущения:

во всех трубах протекает однофазный поток жидкости, температура которого одинакова на всех участках;

все трубы располагаются на одном уровне, в системе нет обратных потоков (рециклов), не учитываются местные сопротивления и перепады давлений в трубах, т.е. рассматриваются так называемые короткие трубопроводы;

системы включают только клапаны (вентили) с постоянными неизменяющимися коэффициентами пропускной способности и закрытые ёмкости (аккумуляторы), давление газа в которых подчиняется идеальным законам.

Рис. 11. Схематическое изображение гидравлической системы с двумя закрытыми емкостями (аккумуляторами) Пример 8. Необходимо смоделировать гидравлическую систему, представленную на рис. 12.

Давления на входах и выходах системы считаются заданными. Также заданы значения коэффициентов пропускной способности. Определить расходы потоков в каждой ветви системы, а также давление в точках разветвления (P11, P12, P13).

Расход через сужающее устройство определяется уравнением:

где P1 – давление до сужающего устройства; P2 – давление после сужающего устройства; k – коэффициент пропускной способности сужающего устройства.

Необходимо добавить 3 уравнения, чтобы система уравнений имела решение (10 переменных – 10 уравнений).

Для выбора алгоритма математической декомпозиции, который позволит определить искомые переменные, необходимо построить и проанализировать информационную матрицу системы уравнений математического описания (48) – (57).

Информационная матрица системы уравнений математического описания представляет собой квадратную матрицу, строки которой соответствуют номерам уравнений, а столбцы – обозначению определяемых переменных.

Составим информационную матрицу (табл. 11).

Каждая строка матрицы соответствует определённому уравнению.

Каждой неизвестной переменной соответствует определённый столбец.

Наличие переменной в том или ином уравнении обозначается символом на пересечении соответствующей строки с соответствующим столбцом.

Данные символы обозначают следующее:

[+] – начальное приближение;

– в данном уравнении переменная определяется при известных значениях остальных переменных;

(+) – переменная косвенно определена другим уравнением.

11. Информационная матрица системы уравнений, описывающей стационарный режим гидравлической системы Если получаемое [U2 (в решении 8) – U2 (в решении 10)] >, то методом половинного деления (или каким-либо другим числовым методом) выбрать новое начальное приближение (для данного примера Р11).

Если U2(8) = U2(10) или U2(8) – U2(10), то модель определена.

Пример 9. Смоделировать барботажный аппарат, схема которого определена на рис. 13.

Условия: давления на входах и выходах известны, также известны коэффициенты пропускной способности.

Определить расходы газа и жидкости через аппарат, давление в аппарате, давление на дне аппарата.

Неизвестные переменные: UG, UL, P, Ph, h,.

где L – плотность жидкости (const); – плотность газожидкостной эмульсии.

Составим информационную матрицу (табл. 12).

12. Информационная матрица системы уравнений, описывающей стационарный режим гидравлической системы UL(3) = UL(4) или UL(3) – UL(4) < – условие окончания алгоритма.

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

1. Определить уравнения расходов через сужающие устройства гидравлической системы, а также давления в точках разветвления системы (варианты контрольных заданий и исходные данные приведены ниже).

2. Закончить полученную систему уравнений для расходов до числа неизвестных давлений в точках разветвления уравнениями для узлов гидравлической системы.

3. Решить полученную систему уравнений на ЭВМ.

ВАРИАНТЫ ЗАДАНИЙ

K1 = 0,75; K2 = 2,2; K3 = 1,9; K4 = 0,3; P1 = 4,7; P2 = 2,8; P3 = 1,9; P4 = 0,5.

K1 = 2,5; K2 = 1,4; K3 = 2; K4 = 1,3; P1 = 7; P2 = 8; P3 = 2,2; P4 = 3,5; P5 = 2.

K11 = 1; K12 = 1,5; K3 = 0,6; K4 = 1,5; K5 = 0,5; K6 = 0,8;

10.

11.

12.

13.

14.

15.

16.

17.

18.

19.

K1 = 0,6; K2 = 1,2; K3 = 2,4; K4 = 1,8; P1 = 4; P2 = 3; P3 = 2; P4 = 1; P5 = 0,5.

20.

K1 = 1,6; K2 = 1,1; K3 = 1,4; K4 = 2,5; K5 = 1; K6 = 1,4; P1 = 5; P2 = 3,6; P3 = 2; P4 = 1,5.

КОНТРОЛЬНЫЕ ВОПРОСЫ

1. Какие основные допущения принимаются при компьютерном моделировании простой гидравлической системы?

2. Как описывается движение потока жидкости через клапан?

3. Как и с какой целью строится информационная матрица системы уравнений математического описания?

4. С какой целью проводится анализ параметрической чувствительности модели?

5. Как провести анализ параметрической чувствительности компьютерной модели?

МОДЕЛИРОВАНИЕ ПРОЦЕССА ПОЛУЧЕНИЯ ПАРА

Цель работы: освоить методику математического моделирования сложных технологических процессов на ЭВМ.

Задания работы:

1. Ознакомиться с основными принципами построения математических моделей на основе балансных уравнений.

2. Выполнить вариант задания.

Содержание отчёта:

1. Вариант задания.

2. Математическая модель процесса.

4. Результаты анализа полученных статических и динамических характеристик объекта.

5. Результаты определения значений настроек ПИ-регулятора АСР.

6. Результаты анализа влияния значений настроек регулятора на качество переходных процессов в АСР.

4. Тексты программ.

КРАТКИЕ ТЕОРЕТИЧЕСКИЕ ПОЛОЖЕНИЯ

Тепловая энергия пара широко используется для отопления, вентиляции и горячего водоснабжения жилых и промышленных помещений.

Для выбора режимов работы, построения системы автоматического регулирования процессами, протекающими при производстве пара, необходимо провести анализ входных и выходных потоков технологического процесса и получить математическое описание объекта управления.

На рисунке 14 представлена упрощённая схема парового котла.

Рис. 14. Упрощённая схема парового котла Проведём анализ процесса производства пара как объекта управления.

Основным качественным выходным параметром процесса производства пара является давление генерируемого пара p.

Входными воздействиями (регулирующими и возмущающими) для данного процесса являются проводимость клапана на линии подачи топлива µ2, проводимость клапана на линии отбора пара µ1 (или расход повх требляемого пара gп), температуры подаваемой топливной смеси t см и питающей воды tв, давление в линии подачи топлива рт и давление в линии потребления пара рпот (рис. 15).

Рис. 15. Связь между входными и выходными параметрами В качестве регулирующего воздействия может быть выбрана проводимость клапана на линии подачи топлива в топку, а остальные воздействия являются возмущающими.

Основным возмущающим воздействием на процесс производства пара является колебание давления в линии потребителя пара.

На основе проведённого анализа разрабатывается математическое описание процесса.

При разработке математического описания приняты следующие допущения.

1. Объём воды и парового пространства считается объектом с сосредоточенными параметрами и принимается идеальное перемешивание в объёме.

2. Тепловой ёмкостью поверхности нагрева пренебрегаем.

3. Паровая среда представляется как идеальный газ.

4. При температуре меньше температуры кипения парообразования не происходит.

5. Потерь тепла через ограждения в окружающую среду не происходит.

6. Температуры воды и пара равны.

7. Удельные теплоёмкости воды, пара, газовоздушной смеси постоянны.

8. Давление в камере сгорания постоянно.

Составим балансные уравнения для процесса производства пара.

Уравнение материального баланса для воды в барабане котла имеет вид:

где Gв – масса воды в барабане котла, кг; g в, g п – расходы входящей питающей воды и испаряющегося пара, соответственно, кг/с.

Уравнение энергетического баланса для воды в барабане котла имеет вид:

где св – удельная теплоемкость воды, Дж/(кг°С); t в – температура воды в барабане котла, °С; tв – температура питающей воды, °С; t к – температура в греющей камере (топке), °С; k – коэффициент теплопередачи через поверхность нагрева, Вт/(м2°С); F – поверхность теплопередачи, м2; g п – расход испаряющегося пара, кг/с; iп – удельная энтальпия испаряющегося пара, Дж/кг.

Продифференцируем (65) как сложную функцию и подставим в выражение (64):

Тогда выражение (65) запишется в виде:

С учётом того, что удельная теплота парообразования равна уравнение энергетического баланса для воды в барабане запишется в виде:

Уравнение материального баланса для пара в котле имеет вид:

где Gп – масса пара в паровом пространстве, кг; g п – расход отбираемого пара, кг/с.

Согласно принятым допущениям пар подчиняется закону идеального газа Менделеева–Клапейрона:

где p – давление пара в котле, Па; Vп – объём парового пространства, м3;

M – молярная масса воды; R – универсальная газовая постоянная.

Объём парового пространства связан с объёмом (или массой воды) в котле соотношением:

где Vo, Vв – общий объём котла и объём воды в котле, м3; в – плотность воды, кг/м3.

Температура кипения воды (и температура пара) является функцией давления:

Расход пара через клапан может быть описан выражением:

где µ1 – проводимость клапана.

Уравнение энергетического баланса для газовоздушной среды в камере сгорания (топке) записывается в виде:

где ск, Gк – удельная теплоёмкость и масса газовоздушной среды в камере сгорания; g см – расход топлива; r – удельная теплота сгорания топлива.

Тепловая емкость камеры сгорания значительно ниже тепловой емкостью нагреваемой воды, поэтому динамикой изменения температуры в топке t к можно пренебречь. В результате уравнение (9) запишется в виде:

Расход топлива также может быть описан зависимостью где µ 2 – проводимость клапана на линии подачи топлива.

Имитационные исследования объекта управления. При автоматизации производственных процессов наиболее правильные решения могут быть приняты на основании данных о статических и динамических свойствах регулируемого объекта. Для получения этих данных определяют соответственно статические и динамические характеристики объекта управления, используя его математическое описание.

Для проведения анализа чувствительности выходных параметров процесса производства пара на изменение входных параметров, а также для выбора регулирующих воздействий проведём построение статических характеристик процесса.

Для этого в математическом описании процесса (64), (66) – (74) приравняем производные в дифференциальных уравнениях к нулю и получим следующую систему алгебраических уравнений:

Для получения статических характеристик объекта управления входные воздействия необходимо менять в диапазоне ±20% от своего значения в номинальном режиме.

Графики статических характеристик приведены на рис. 16.

Из анализа статических характеристик можно сделать вывод, что давление обладает наибольшей чувствительностью по изменению проводимости клапанов на трубопроводах подвода топливной смеси и отбора пара из парового котла. Следовательно, проводимость клапанов можно использовать в качестве регулирующего воздействия. Давление в линии подачи топливной смеси можно считать постоянным. Основным возмущением является давление потребляемого пара, так как чувствительность давления пара в котле на изменение этого параметра наибольшая.

Анализ динамических характеристик котла проведём по переходным функциям, построенным для наиболее сильнодействующего возмущения pпот и регулирующего воздействия µ 2 (рис. 17). При этом значения входных воздействий менялись на 5% от своего значения в номинальном режиме.

Рис. 16. Статические характеристики парового котла:

Рис. 17. Переходные функции парового котла при изменении:

Моделирование одноконтурной системы регулирования давления в паровом котле. Анализируя переходные функции можно сделать следующие выводы:

объект обладает свойством самовыравнивания;

время перехода из одного состояния в другое составляет 1300 с.

Структурная схема одноконтурной АСР давления пара представлена на рис. 18.

Рис. 18. Структурная схема системы регулирования давления пара:

х – вход объекта; p – давление – выходная регулируемая величина;

pз – заданное значение температуры; pи – измеренное значение температуры;

Помимо собственно объекта регулирования в систему регулирования входят чувствительный элемент (Д), исполнительный механизм с регулирующим органом (ИМ) и регулятор (R). Поэтому в математическое описание АСР входят уравнения, описывающие характеристики датчика, исполнительного механизма и регулятора.

Динамические характеристики датчика и исполнительного механизма представляются в виде дифференциальных уравнений:

где TД, TИМ – постоянные времени датчика давления и исполнительного мость клапана на линии подачи топлива.

Постоянные времени используемого датчика давления ТД = 0,5 мин, а исполнительного механизма ТИМ = 1 мин.

Расчёт одноконтурной АСР заключается в определении настроечных параметров функции R1 ( p ). В практике автоматизации промышленных объектов широкое распространение получили ПИ-регуляторы. Для регулирования давления в паровом котле уравнение ПИ-регулятора имеет вид:

где K р, Tи – параметры настройки регулятора; е – ошибка регулирования.

Ошибка регулирования определяется как разница между заданным значением и текущим значением измеренного давления в паровом котле:

При заданном законе регулирования расчёт одноконтурной АСР сводится к определению оптимальных (с точки зрения какого-либо критерия качества) настроечных параметров регулятора с учётом выполнения ограничения на запас устойчивости системы.

Наиболее часто в качестве критерия качества используется интегральный квадратичный критерий при ограничении на степень затухания:

где кон – время регулирования.

На рисунках 19, 20 показаны графики переходных процессов в одноконтурной АСР при найденных настройках регулятора. При этом задание регулятору давления пара в парогенераторе изменялось с 8,995 105 Па до 9,995 105 Па. По каналу возмущения – давление пара в линии потребления изменялось с 4,905 105 Па до 6,867 105 Па.

Рис. 19. Переходной процесс в АСР при изменении задания Рис. 20. Переходной процесс в АСР при действии возмущения Анализ качественных показатели переходных процессов в АСР.

Качество регулирования численно может быть охарактеризовано прямыми показателями качества, которые определяются непосредственно из графиков переходных процессов АСР.

В табл. 13 сведены прямые показатели качества процесса регулирования в АСР (с ПИ-регулятором).

13. Прямые показатели качества процесса регулирования Программа построения переходных процессов в АСР давления пара в котле представлена ниже.

ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

1. В соответствии с заданным вариантом исходных данных процесса получения пара в паровом котле (варианты заданий выдаются преподавателем индивидуально) необходимо исследовать статические и динамические характеристики процесса.

2. Определить значения параметров регулятора и произвести построение переходных процессов в АСР.

3. Исследовать влияние значения параметров настройки регулятора на качество переходных процессов. Для этого значения настроек регулятора необходимо менять относительно найденных значений в большую и меньшую сторону.

Программа для расчёта уравнений математической модели и построения переходных процессов в АСР dt=0.05;

tk=1000;

t=[0:dt:tk];

n=length(t);

t_v=zeros(1,n);

t_k=zeros(1,n);

p=zeros(1,n);

M_v=zeros(1,n);

M_p=zeros(1,n);

G_p=zeros(1,n);

pi=zeros(1,n);

% начальные условия t_v(1)=173.9754;

t_k(1)=278.1349;

p(1)=8.9951e+005;

M_v(1)=2.5e+003;

pi(1)=p(1);

M_p(1)=32.6929;

% описание теплофизических характеристик c_v=1*4.19e3;

c_k=1500;

M=18e-3;

R=8.31;

ro_p=8.4;

ro_k=1;

ro_v=1000;

r_t=3.6e7;

Lamda=1945e3;

% описание характеристик парового котла и значений входов V_0=10;

K=1500;

F=10;

G_v_n=0.6032;

t_v_vx=20;

t_t_vx=20;

p_t=2*9.81e4;

p_pot=5*9.81e4;

p_k=1*9.81e4;

% проводимости клапанов mu1=0.0000008;

mu2_n=0.0000002;

% характеристики насыщенного пара xx=[0.98 0.99 1 1.01 1.013 1.02 2 3 4 5 6 7 8 9 10 20 40]*1e5;

yy=[99.23 99.46 99.7 99.93 100 100.27 120 133 143 151 158 164 169 212 250];

% переменные для определения интегральных составляющих ПИ-регуляторов sum=0;

sum2=0;

% настройки ПИ-регуляторов Kp2=0.000000000002;

Ti2=100;

Kp=20;

Ti=20;

p_z=9.9951e+005;

M_v_z=2500;

% характеристики датчика и ИМ Td=10;

Tim=20;

% цикл решения уравнений модели методом Эйлера for i=2:n % поддержание номинальной массы воды в котле eps=M_v(i-1)-M_v_z; % ошибка регулирования sum=sum+eps;

G_v=G_v_n-Kp*(eps+sum*dt/Ti); % уравнение регулятора if G_v=5;

G_v=5;

% определение давление пара по температуре воды p(i-1)=interp1(yy,xx,t_v(i-1),'linear');

% определение измеренного значения давления pi(i)=pi(i-1)+dt*((p(i-1)-pi(i-1))/Td);

eps2=-pi(i)+p_z; % ошибка регулирования sum2=sum2+eps2;

mu2_=mu2_n+Kp2*(eps2+sum2*dt/Ti2); % уравнение регулятора % ограничение на проводимость клапана if mu2_ < mu2_=0;

if i== mu2_old=mu2_;

% определение реальной проводимости клапана mu2=mu2_old+dt*(mu2_-mu2_old)/Tim;

mu2_old=mu2;

% решение уравнений модели V_p=V_0-M_v(i-1)/ro_v;

M_p(i)=p(i-1)*V_p*M/(R*(t_v(i-1)+273));

if p_pot > p(i-1) G_p(i-1)=0;

else G_p(i-1)=mu1*sqrt(p(i-1)^2-p_pot^2);

end G_pi=G_p(i-1)+(M_p(i)-M_p(i-1))/dt;