«МОДЕЛИРОВАНИЕ СИСТЕМ Издательство ФГБОУ ВПО ТГТУ Учебное издание ЕЛИЗАРОВ Игорь Александрович, МАРТЕМЬЯНОВ Юрий Федорович, СХИРТЛАДЗЕ Александр Георгиевич, ТРЕТЬЯКОВ Александр Александрович МОДЕЛИРОВАНИЕ СИСТЕМ Учебное ...»
И.А. Елизаров
Ю.Ф. Мартемьянов
А.Г. Схиртладзе
А.А. Третьяков
МОДЕЛИРОВАНИЕ
СИСТЕМ
Издательство ФГБОУ ВПО «ТГТУ»
Учебное издание
ЕЛИЗАРОВ Игорь Александрович,
МАРТЕМЬЯНОВ Юрий Федорович,
СХИРТЛАДЗЕ Александр Георгиевич,
ТРЕТЬЯКОВ Александр Александрович
МОДЕЛИРОВАНИЕ СИСТЕМ
Учебное пособие Редактор Т.М. Г л и н к и н а Инженер по компьютерному макетированию М.А. Ф и л а т о в а Подписано в печать 19.09.2011 Формат 60 84/16. 5,58 усл. печ. л. Тираж 100 экз. Заказ № 384 Издательско-полиграфический центр ФГБОУ ВПО «ТГТУ»392000, г. Тамбов, ул. Советская, д. 106, к. Министерство образования и науки Российской Федерации Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования «Тамбовский государственный технический университет»
И.А. ЕЛИЗАРОВ, Ю.Ф. МАРТЕМЬЯНОВ, А.Г. СХИРТЛАДЗЕ, А.А. ТРЕТЬЯКОВ
МОДЕЛИРОВАНИЕ СИСТЕМ
Допущено Учебно-методическим объединением вузов по образованию в области автоматизированного машиностроения в качестве учебного пособия для студентов высших учебных заведений, обучающихся по направлению подготовки «Автоматизация технологических процессов и производств»Тамбов Издательство ФГБОУ ВПО «ТГТУ»
УДК 519. ББК 22. Е Рецензенты:
Кандидат технических наук, доцент ФГБОУ ВПО «ТГТУ»
О.Г. Иванова Кандидат технических наук, доцент, первый заместитель генерального директора ОАО «Корпорация «Росхимзащита»
С.Б. Путин Е511 Моделирование систем : учебное пособие / И.А. Елизаров, Ю.Ф. Мартемьянов, А.Г. Схиртладзе, А.А. Третьяков. – Тамбов :
Изд-во ФГБОУ ВПО «ТГТУ», 2011. – 96 с. – 100 экз.
ISBN 978-5-8265-1022- Содержит серию лабораторных работ для практического освоения методов математического моделирования и приобретения навыков проведения вычислительных (компьютерных) расчётов при создании систем и средств моделирования и соответствующих математических расчётов, работы по построению математических моделей на основе данных активного эксперимента, исследованию линейных стационарных систем управления, построению аналитических и экспериментально-аналитических моделей ряда типовых технологических процессов. Сопровождается необходимыми теоретическими положениями для выполнения лабораторных работ с использованием системы автоматизированных расчётов MATLAB+SIMULINK.
Предназначено для студентов высших учебных заведений, обучающихся по направлениям подготовки «Автоматизация технологических процессов и производств».
УДК 519. ББК 22. © Федеральное государственное бюджетное ISBN 978-5-8265-1022- образовательное учреждение
ПРЕДИСЛОВИЕ
В подготовке студентов, специализирующихся в области построения автоматизированных систем управления технологическими процессами, одно из ведущих мест занимает курс «Моделирование систем». Эффективность освоения этого курса в значительной мере зависит от содержания и постановки лабораторного практикума.В общем случае под термином «математическое моделирование автоматических систем» понимаются процессы отыскания их математических моделей, а также непосредственного исследования и анализа этих моделей на основе методов теории автоматического управления аналитически, графоаналитически или с использованием компьютеров. По мере развития вычислительной техники разработано достаточно много прикладных программ, обеспечивающих исследование переходных и установившихся процессов в автоматических системах любой сложности и практически в любых режимах работы при изменении их параметров и структуры. Среди них – пакеты MATLAB/SIMULINK, MATHCAD, VISSIM, LABVIEW и др.
Система MATLAB разработана фирмой The MathWorks, Inc. (США) и является системой инженерных и научных расчётов для различных областей науки и техники. Среди них:
математика и вычисления;
разработка алгоритмов;
вычислительный эксперимент, имитационное моделирование;
анализ данных, исследование и визуализация результатов;
исследования в области автоматического управления;
статистическая обработка сигналов и процессов и др.
1. Интерактивная система, которая позволяет производить вычисления с помощью непосредственного ввода команд с клавиатуры;
2. Огромная библиотека готовых функций и алгоритмов, реализующих наиболее распространённые методы вычислений, которые для облегчения их поиска разбиты на специализированные разделы (Toolboxes);
3. Язык программирования, позволяющий пользователю создавать собственные функции, а также законченные программные приложения, в том числе и с использованием библиотечных функций системы.
Программа SIMULINK является расширением пакета MATLAB, достаточно самостоятельным его инструментом. При работе с SIMULINK практически не требуется знать сам MATLAB и остальные его приложения.
С другой стороны, доступ к функциям MATLAB и другим его инструментам остаётся открытым и их можно использовать в SIMULINK.
Часть входящих в состав пакетов имеет инструменты, встраиваемые в SIMULINK (например, LTI-Viewer приложения Control System Toolbox – пакета для разработки систем управления).
Как программное средство SIMULINK – типичный представитель визуально-ориентированных языков программирования. На всех этапах работы, особенно при подготовке моделей систем, пользователь практически не имеет дела с обычным программированием. Программа в кодах автоматически генерируется в процессе ввода выбранных блоков компонентов, их соединений и задания параметров компонентов.
При работе с SIMULINK пользователь имеет возможность модернизировать библиотечные блоки, создавать свои собственные, а также составлять новые библиотеки блоков. При моделировании пользователь может выбирать метод решения дифференциальных уравнений. В ходе моделирования имеется возможность следить за процессами, происходящими в системе. Для этого используются специальные устройства наблюдения, входящие в состав библиотеки SIMULINK. Результаты моделирования могут быть представлены в виде графиков или таблиц.
1. КРАТКОЕ ОПИСАНИЕ РАБОТЫ В MATLAB
Работа MATLAB возможна в двух режимах: командном режиме и режиме выполнения m-файла.Работа системы в командном режиме напоминает работу в качестве калькулятора. В командной строке MATLAB после символа » можно сразу вводить исходные данные для вычисления с помощью встроенного текстового редактора. Ввод данных завершается нажатием клавиши ENTER.
Результаты вычислений выводятся с новой строки (без маркера »).
Для блокировки вывода результатов вычислений некоторого выражения его надо закончить знаком ; (точка с запятой). Это бывает удобным для скрытия результатов промежуточных вычислений.
В качестве знака присваивания используется знак равенства =. Например, ввод строки x = 1 означает объявление и инициализацию переменной x.
Имя переменной может содержать до 30 символов и должно не совпадать с именами функций, процедур системы и системных переменных.
При этом система различает большие и маленькие буквы в переменных, т.е. переменные ааа и ААА в MATLAB обозначают разные переменные.
Система MATLAB имеет несколько имен переменных, которые используются самой системой и входят в состав зарезервированных:
i, j – мнимая единица (корень квадратный из –1);
pi – число (сохраняется в виде 3.141592653589793);
inf – обозначение машинной бесконечности;
NaN – обозначение неопределённого результата (например, типа 0/0 или inf/inf);
eps – погрешность операций над числами с плавающей запятой;
ans – результат последней операции без знака присваивания;
realmax и realmin – максимально и минимально возможные величины числа, которые могут быть использованы.
Эти переменные можно использовать в математических выражениях.
Встроенные функции записываются строчными буквами, а их аргументы приводятся в круглых скобках, например, sin(x).
В одном сеансе работы системы MATLAB можно определить несколько переменных и/или вычислить несколько выражений, например:
Для вывода значения конкретной переменной достаточно ввести её имя, например:
Если вывод результата вычисления только что введённого выражения не заблокирован и не указана переменная, значение которой необходимо вывести, то система MATLAB сама назначает такую переменную с именем ans, присваивает ей значение последнего выражения и выводит её как результат вычислений, например:
» sin(pi/2) Система MATLAB изначально создавалась как «матричный вычислитель» и оптимизирована для проведения вычислений с векторами, матрицами и массивами. Более того, в системе по умолчанию предполагается, что каждая переменная – это вектор или матрица. Например, если задано X = 1, то система воспринимает это как задание вектора с единственным элементом, значение которого равно 1. В свою очередь, вектор – это матрица, число столбцов которой равно 1.
Для задания вектора с большим числом элементов, их значения надо перечислить в квадратных скобках, используя в качестве разделителя пробел или запятую. Задание матриц осуществляется аналогично, для разделения строк используется знак ; (точка с запятой) Например:
В MATLAB предусмотрено наличие большого количества операций, которые можно совершать с матрицами и векторами. Среди них:
• сложение, вычитание (+,–);
• умножение (*);
• обращение (inv);
• деление слева направо (/) и справа налево (\);
• возведение в степень (^);
• транспонирование (‘);
• создание нижней треугольной матрицы А: tril(A);
• создание верхней треугольной матрицы А: triu(A);
• формирование единичной матрицы заданного размера n: eye(n);
• формирование матрицы единиц заданного размера n m :
ones(n,m). Для создания квадратной матрицы: ones(n);
• формирование матрицы нулей заданного размера n m : zeros(n,m). Для создания квадратной матрицы: zeros(n);
• извлечение диагонали заданной матрицы А: daig(A);
• вычисление следа матрицы А: trace(A);
• собственные числа матрицы А: eig(A);
В MATLAB предусмотрено наличие большого числа встроенных математических функций [1]. Среди них:
– элементарные функции:
• sin(Z), cos(Z), tan(Z), cot(Z), sec(Z), csc(Z) – синус, косинус, тангенс, котангенс числа Z;
• sinh(Z), cosh(Z), tanh(Z), coth(Z), sech(Z), csch(Z) – гиперболические синус, косинус, тангенс, котангенс, секанс, косеканс числа Z;
• asin(Z), acos(Z), atan(Z), acot(Z), asec(Z), acsc(Z) – арксинус (в радианах в диапазоне от -/2 до +/2), арккосинус (в диапазоне от 0 до ), арктангенс (в диапазоне от -/2 до +/2), арккотангенс, арксеканс, арккосеканс;
• asinh(Z), acosh(Z), atanh(Z), acoth(Z), asech(Z), acsch(Z) – обратные гиперболические синус, косинус, тангенс, котангенс, секанс, косеканс;
• exp(Z) – экспонента числа Z;
• log(Z) – натуральный логарифм;
• log10(Z) – десятичный логарифм;
• sqrt(Z) – квадратный корень из числа Z;
• abs(Z) – модуль числа Z.
• round(Z) – обычное округление числа Z к ближайшему целому;
• mod(X,Y) – целочисленное деление X на Y;
• rem(X,Y) – вычисление остатка от деления X на Y;
• sign(Z) – вычисление сигнум-функции числа Z;
– функции комплексного аргумента.
Практически все вышеперечисленные элементарные математические функции вычисляются при комплексных значениях аргумента и получают в результате этого комплексные значения результата. Например, функция sqrt вычисляет, в отличие от других языков программирования, квадратный корень из отрицательного аргумента, а функция abs при комплексном значении аргумента вычисляет модуль комплексного числа.
В MATLAB есть несколько дополнительных функций, рассчитанных только на комплексный аргумент:
• real(Z) – выделяет действительную часть комплексного аргумента Z;
• imag(Z) – выделяет мнимую часть комплексного аргумента;
• angle(Z) – вычисляет значение аргумента комплексного числа Z (в радианах в диапазоне от – до +);
• conj(Z) – выдает число, комплексно сопряжённое относительно Z.
1.1.3. Решение системы линейных алгебраических уравнений Рассмотрим систему линейных алгебраических уравнений:
В матричной форме эта система может быть записана в виде AX = B, где матрицы A, B, X имеют вид:
Решение системы уравнений в матричной форме записи представляется выражением:
где A – матрица, обратная матрице A.
В MATLAB это решение может быть найдено при выполнении выражений:
» X=inv(A)*B или П р и м е р. С использованием MATLAB найти решение системы линейных уравнений:
Решение:
Оператор : (двоеточие) используется для создания упорядоченной последовательности чисел с равноотстоящими значениями. Формат записи этого оператора:
[Начальное_значение]:[ Шаг]:[Конечное_значение] Если шаг не задан, то он принимается равным 1 или –1 в зависимости от того, больше конечное значение начального значения или меньше.
При его использовании создаётся вектор:
» k = 2:2: » s=sin(m) При необходимости осуществления почленного умножения (или деления) элементов одного массива на элементы другого массива (той же размерности) используются соответственно операторы (.*) и (./). Например:
Для осуществления операции возведения в степень элементов массива используется операция ".^ ".
Одной из самых привлекательных черт системы MATLAB является визуализация вычислений. Для построения графиков функций может использоваться функция plot.
Эта функция имеет различные формы записи в зависимости от набора входных аргументов. Если y – вектор, функция plot(y) производит построение кусочно-линейного графика элементов y в зависимости от индекса элементов y. Если функция plot(x,y) имеет два векторных аргумента одинаковой размерности, то будет построен график y в зависимости от x.
Функция plot может иметь несколько пар аргументов x – y. В этом случае она изображает несколько функций y(x) кривыми разного цвета, используя предопределённый (но изменяемый пользователем) список цветов. Например, plot(x1, y1, x2, y2, x3, y3) – выстраиваются три графика функции: y1(x1), y2(x2), y3(x3).
Также имеется возможность определять цвет, стиль линии и маркеры с помощью ещё одного аргумента функции plot:
где color_style_marker – это символьная строка (ограниченная одиночными кавычками), характеризующая цвет, стиль линии и тип маркера:
• символы определения цвета: 'c', 'm', 'y', 'r', 'g', 'b', 'w', 'k'. Они соответствуют цветам cyan, magenta, yellow, red, green, blue, white, black;
• обозначения стиля линии: '-' сплошная, '--' пунктирная, ':' штриховая, '-.' штрихпунктирная;
• типы маркера: '+', 'o', '*', 'x'.
Например, выражение: plot(x, y, 'r--o') вычерчивает красную пунктирную кривую и устанавливает маркеры «кружок» в каждой точке данных.
Любая из трёх перечисленных составляющих color_style_marker может быть пропущена. Если определяется тип маркера, а стиль линии опущен, система MATLAB строит только маркеры для этой кривой.
На график дополнительно можно нанести сетку из координатных линий, информацию о кривой (заголовок графика), а также обозначить координатные оси.
Нанесение координатной сетки на график осуществляется с помощью функции grid, к которой следует обратиться сразу после обращения к функции plot:
» plot(x, y), grid Заголовок графика выводится с помощью процедуры title. Если после обращения к процедуре plot вызвать title таким образом:
то над графиком появится текст, записанный между апострофами в скобках. При этом следует помнить, что текст всегда должен помещаться в апострофы.
Аналогично можно вывести объяснения к графику, которые размещаются вдоль горизонтальной оси (функция xlabel) и вдоль вертикальной оси (функция ylabel).
Например, совокупность операторов » x = -4*pi : pi/100 : 4*pi;
» y = 2* sin(x+pi/3);
» plot(x,y), grid;
» title('Функция y = 2*sin(x+pi/3)');
» xlabel('x'); ylabel('y');
приведёт к оформлению поля фигуры в виде:
1.2. ПРОГРАММИРОВАНИЕ В СРЕДЕ MATLAB
Работа в режиме калькулятора в среде MATLAB, несмотря на довольно значительные возможности, во многих отношениях неудобна. Невозможно повторить предшествующие вычисления и действия при новых значениях исходных данных без повторного набора предшествующих операторов. Нельзя возвратиться назад и повторить некоторые действия, или по некоторому условию перейти к выполнению другой последовательности операторов.Поэтому сложные алгоритмы вычислений, сопровождаемые переходами по определённым условиям с часто повторяемыми однотипными действиями, рациональнее оформлять в виде программ.
Создание программы в среде MATLAB осуществляется с помощью либо собственного встроенного, либо стороннего текстового редактора (например, Notepad среды Windows).
Программы на языке MATLAB имеют две разновидности – так называемые Script-файлы (файлы-сценарии, или управляющие программы) и файлы-функции (процедуры). Обе разновидности должны иметь расширение имени файла.m. С помощью Script-файлов оформляют основные программы, управляющие от начала до конца организацией всего вычислительного процесса, и отдельные части основных программ (они могут быть записаны в виде отдельных Script-файлов). Как файл-функции оформляются отдельные процедуры и функции (т.е. такие части программы, которые рассчитаны на неоднократное использование Script-файлами или другими процедурами при измененных значениях исходных параметров и не могут быть выполнены без предварительного задания значений переменных, которые называют входными) [1].
В дальнейшем под М-файлом будем понимать любой файл (файлфункцию или Script-файл), записанный на языке системы MATLAB.
Основные особенности записи текста программы (М-файла) [1]:
Обычно каждый оператор записывается в отдельной строке текста программы. Окончание ввода строки осуществляется нажатием клавиши.
Можно размещать несколько операторов в одной строке, которые отделяются друг от друга символом ' ; ' или ', '. Можно длинный оператор записывать в несколько строк. При этом предыдущая строка оператора должна заканчиваться тремя точками '... '.
Если очередной оператор не заканчивается символом ' ; ', результат его действия при выполнении программы будет выведен в командное окно. Чтобы предотвратить вывод на экран результатов действия оператора программы, запись этого оператора в тексте программы должна заканчиваться символом ' ; '. Если выражение завершается знаком ' ; ', то MATLAB не выводит результирующее значение на экран. Обычно формат оператора с ' ; ' используется при работе со средой в интерактивном режиме, а без ' ; ' – при написании программы.
Строка программы, начинающаяся с символа ' % ', не выполняется. Эта строка воспринимается системой MATLAB как комментарий.
Строки комментария, предшествующие первому выполняемому оператору программы, т.е. такому, который не является комментарием, воспринимаются системой MATLAB как описание программы.
В программах на языке MATLAB отсутствует символ окончания текста программы.
В языке MATLAB переменные не описываются и не объявляются.
Любое новое имя, появляющееся в тексте программы при ее выполнении, воспринимается системой MATLAB как имя матрицы. Размер этой матрицы устанавливается при предварительном вводе значений ее элементов либо определяется действиями по установлению значений ее элементов, описанными в предшествующих операторах или процедуре.
В языке MATLAB невозможно использование матрицы или переменной, в которой предварительно не введены или не вычислены значения ее элементов (при выполнении программы MATLAB появится сообщение об ошибке – "Переменная не определена").
1.2.1. Создание простейших файл-функций Файл-функция (процедура) должна начинаться со строки заголовка:
где, – соответственно, перечни выходных и входных величин файл-функции. Если перечень выходных величин (OutVar) содержит только один объект (в общем случае – матрицу), то файл-функция представляет собой обычную функцию (одной или нескольких переменных). Первая строка в этом случае имеет вид:
function = (< InVar >).
Если же в результате выполнения файл-функции должны быть определены (вычислены) несколько объектов (матриц), такая файл-функция представляет собой уже более сложный объект. Общий вид первой строки в этом случае становится таким:
function [y1, y2,..., y] = (< InVar>), т.е. перечень выходных величин y1, y2,..., y должен быть представлен как вектор-строка с элементами y1, y2,..., y (все они могут быть матрицами).
В простейшем случае функции одной переменной заголовок приобретет вид:
function y = func(x), где func – имя функции (М-файла).
Пример: составить m-файл для вычисления функции Процесс создания m-файла выглядит следующим образом:
1. Вызвать меню File командного окна MATLAB и выбрать в нём сначала команду New, а затем команду M-file.
2. В появившемся окне текстового редактора набрать текст:
function y = func(x,z) y=log(x+z)+x./(x+z);
3. Сохранить этот текст в файле под именем func.m. Необходимый М-файл создан.
Теперь можно пользоваться этой функцией при расчётах. Например:
При создании файл-функций необходимо помнить некоторые особенности:
• Во избежание вывода на экран промежуточных результатов, необходимо в тексте процедуры все вычислительные операторы завершать символом " ; ".
• С точки зрения области видимости переменных переменные в файл-функции являются локальными, поэтому имена переменных, используемые в файл-функции, могут не совпадать с именами соответствующих переменных при обращении к этой файл-функции.
• С целью использования в файл-функции глобальной переменной (некоторой переменной рабочего пространства (Workcpace)) необходимо в диалоговом режиме (или в Sсript-файле) и файл-функции объявить эту переменную с атрибутом global. Если в одной строке объявляются несколько переменных как глобальные, они должны отделяться пробелами.
• Возможность использования файл-функции как для отдельных чисел, так и для векторов и матриц обусловлена применением в записи соответствующего М-файла вместо обычных знаков арифметических действий их аналогов с предшествующей точкой '. ' ('.* ', './ ', '.^ ' и др.).
При создании Script-файлов учитываются следующие особенности [1]:
• Script-файлы являются независимо (самостоятельно) исполняемыми блоками операторов и команд;
• все используемые переменные образуют так называемое рабочее пространство (Workspace), которое является общим для всех исполняемых Script-файлов; из этого следует, что при выполнении нескольких Scriptфайлов имена переменных в них должны быть согласованы, так как одно имя означает в каждом из них один и тот же объект вычислений;
• в них отсутствует заголовок, т.е. первая строка определённого вида и назначения (в отличие от файл-функций, которые начинаются с заголовка function);
• обращение к ним не требует указания никаких имен переменных:
все переменные формируются в результате выполнения программы либо сформированы ранее и существуют в рабочем пространстве.
1.2.3. Ввод и вывод информации в диалоговом режиме Для обеспечения взаимодействия с пользователем в процессе выполнения m-файла используются такие команды:
disp, sprintf, fprintf, input, menu, keyboard, pause.
Команда disp осуществляет вывод значений указанной переменной или указанного текста в командное окно. Обращение к ней имеет вид:
disp(x1) disp(‘value’) Чтобы вывести значения нескольких переменных в одну строку (например, при создании таблиц данных), нужно создать единый объект, который содержал бы все эти значения. Это можно сделать, объединив соответствующие переменные в вектор, пользуясь операцией создания векторастроки:
Тогда вывод значений нескольких переменных в одну строку будет иметь вид:
Для одновременного вывода символьной и цифровой информации в командное окно удобно использовать функцию sprintf, например:
» disp(sprintf('Параметр1 = %g Параметр2= %g',x1,x2)) Параметр1 = –3.14 параметр2= –2. Аналогично работает функция fprintf, где в качестве первого параметра для вывода информации в командное окно используется 1, например:
» fprintf(1,'Параметр1 = %g Параметр2= %g',x1,x2) Ввод информации с клавиатуры в диалоговом режиме можно осуществить с помощью функции input:
При выполнении этой функции программа ожидает ввода информации с клавиатуры. По окончании ввода, которое определяется нажатием клавиши, введённая информация запоминается в программе под именем "х", и выполнение программы продолжается.
Команда pause временно прекращает выполнение программы до тех пор, пока пользователь не нажмёт любую клавишу клавиатуры.
В языке MATLAB для организации ветвления используются команды if и switch.
Общая форма записи оператора if имеет вид [elseif ] else end Пример. Разработать m-файл-функцию, осуществляющую определение знака числа.
Функция будет располагаться в m-файле signum.m.
function у = signum(x) if х > elseif х == else end Команда switch позволяет осуществлять ветвление по нескольким условиям (направлениям) так же просто, как и с двумя, причём условия рассматриваются на равенство. Ниже представлен простой пример с разделением на три условия формирования параметра ввода.
function у = count(х) switch х case case otherwise end Здесь выражение switch вычисляет параметр ввода х, а затем выполнение программы в файле переходит туда, где выражение case имеет то же самое значение, что и параметр ввода. Таким образом, если параметр ввода х имеет значение, равное 1, тогда параметр вывода у определяется в виде строки 'one', если х равен 2, тогда у определяется в виде строки 'two'.
После выполнения всех команд, следующих за выбранным условием, программа MATLAB встречает либо другое выражение case, либо выражение otherwise, что приводит к переходу выполнения программы к выражению end.
В языке MATLAB для организации циклов используются команды for и while.
Организация цикла с помощью команды for выглядит следующим образом:
for n=n_0: [step:] n_k end Цикл начинается с выражения for и заканчивается выражением end.
Величины n_0 – начальное значение переменной цикла, step – шаг, n_k – конечное значение переменной цикла.
Пример. Требуется вычислить 10!.
Фрагмент кода для определения факториала представлен ниже. Он может быть выполнен непосредственно в командном режиме f=1;
for n=2:1: end Организация цикла с помощью команды while выглядит следующим образом:
while end Цикл while продолжается до тех пор, пока истинно.
Пример: вычислить ряд n = 1;
oldsum = -1;
newsum = 0;
while (newsum-oldsum)>1e- n = n +1;
end newsum newsum = 1. Фрагмент кода для вычисления ряда представлен ниже. Он может быть выполнен непосредственно в командном режиме.
Иногда бывает необходимо, чтобы программа MATLAB преждевременно вышла из цикла, например, при возникновении определённого условия. Для этой цели используется команда break.
2. КРАТКОЕ ОПИСАНИЕ РАБОТЫ В SIMULINK
SIMULINK – это интерактивная система для анализа линейных и нелинейных динамических систем. Она может работать с линейными, нелинейными, непрерывными, дискретными, многомерными системами.Перед построением модели необходимо предварительно загрузить систему MATLAB и запустить пакет SIMULINK.
Процесс построения модели SIMULINK включает в себя компоновку модели и задание необходимых параметров. Компоновка заключается в выборе из библиотек SIMULINK необходимых блоков, их размещение в открывшемся окне и соединение между собой. Далее для каждого блока устанавливаются соответствующие параметры, отвечающие требованиям моделируемой системы. Для их изменения надо дважды щелкнуть на блоке и изменить нужные значения в диалоговом окне.
Блоки имеют названия, которые можно изменять, щёлкнув по нему левой клавишей мыши и отредактировав текст.
Блоки можно поворачивать. Для того, чтобы повернуть блок на 90 градусов, надо выделить его и нажать клавиши Ctrl+R. Комбинация Ctrl+I позволяет выполнить зеркальное отражение входов и выходов.
Блоки соединяются линиями связи, по которым распространяются сигналы. Для того, чтобы соединить блоки, надо щёлкнуть левой клавишей мыши по источнику сигнала и затем, при нажатой клавише Ctrl, по блоку-приёмнику. Можно также протянуть мышкой линию связи между нужными выходом и входом [2].
При создании ответвления от существующей соединительной линии, т.е. для соединения входного порта какого-либо блока с существующей линией, необходимо при выполнении соединения воспользоваться клавишей Ctrl. Также для создания точки разветвления в соединительной линии можно подвести курсор к предполагаемому узлу и, нажав правую клавишу мыши, протянуть линию.
Модель можно скопировать в буфер обмена в виде растрового рисунка. Для этого в окне модели надо выбрать в верхнем меню пункт Edit / Copy model to clipboard. Предварительно лучше уменьшить размеры окна до минимальных, чтобы не было белых полей.
Блок-диаграммы SIMULINK могут быть объединены в составные блоки, что позволяет использовать иерархическое представление структуры модели. Создание составного блока возможно двумя путями:
1. Перенести из библиотеки блоков блок Subsystem (составной блок) и, раскрыв его окно, сформировать его структуру в этом отдельном окне.
2. Выстроить блок-диаграмму в основном окне создаваемой модели, выделить её часть (предназначенную для объединения в составной блок) и, перейдя в меню, выбрать пункт Edit/Create Subsystem. В результате в основном окне автоматически образуется составной блок, в котором будет заключена та часть блок-диаграммы, что была выделена.
Для того, чтобы запустить моделирование, надо щёлкнуть левую клавишу мыши по кнопке на панели инструментов. Эта же кнопка позволяет остановить моделирование при необходимости. Запуск модели на расчёт можно осуществить также через пункт меню Simulation/Start.
Предварительно перед запуском разработанной модели необходимо:
1. Установить параметры работы модели (интервал времени работы модели, метод интегрирования, погрешность и др.). Установка этих параметров осуществляется через пункт меню Simulation/Parameters… Самые важные параметры – это время моделирования (Stop time) и метод численного интегрирования уравнений (Solver options).
2. Сохранить разработанную модель (пункт File/Save).
Библиотеки блоков SIMULINK. Доступные блоки SIMULINK отображаются в окне Simulink Library Browser. Все блоки по функциональному назначению разделены по библиотекам. Среди них:
• Continuous – библиотека непрерывных элементов (интегратор, дифференциатор, линейная система обыкновенных дифференциальных уравнений и др.);
• Discontinuities – нелинейные элементы (реле, зона нечувствительности, люфт и др.);
• Discrete – библиотека дискретных элементов (интегратор с дискретным временем, дискретный фильтр и др.);
• Math operations – математические функции;
• Logic and Bit operations – логические и битовые операции;
• Sinks – средства отображения (временная диаграмма, вывод результатов в файл, остановка выполнения модели и т.д.);
• Sources – источники сигналов (генератор импульсных/синусоидальных сигналов, генератор случайных чисел, генератор пилообразных сигналов, часы и т.д.).
Непрерывные линейные системы (Continuous) Transfer Fcn – передаточная функция, в параметрах задаются числитель (Numerator) и знаменатель (Denominator) в виде полиномов.
State Space – модель в пространстве состояний, в параметрах задаётся четвёрка матриц, определяющих модель, и начальные условия для вектора состояния Zero-Pole – модель в форме «нули-полюса», в параметрах задаются массивы нулей (Zeros), полюсов (Poles), а также коэффициент усиления (Gain).
Integrator – интегратор с возможностью установки начальных условий (Initial condition), а также пределов насыщения (Lower saturation limit и Upper saturation limit). Когда сигнал выхода выходит за границы, определяемые этими пределами, интегрирование Transport Delay – блок фиксированной задержки сигнала.
Нелинейные элементы (Discontinuities) Relay – блок реализует релейную характеристику с Saturation – блок выполняет ограничение величины Backlash – блок моделирует нелинейность типа Rate Limiter – блок обеспечивает ограничение скорости изменения сигнала (первой производной).
Dead Zone – блок реализует нелинейную зависимость типа «зона нечувствительности (мертвая зона)».
Основные источники сигналов (Sources) Constant – сигнал постоянной величины.
Step – ступенчатый сигнал, меняется время скачка (Step Time), начальное (Initial Value) и конечное значения (Final Ramp – линейно возрастающий сигнал с заданным наклоном (Slope). Можно задать также время начала изменения сигнала (Start Time) и начальное значение (Initial Value).
Pulse Generator – генератор прямоугольных импульсов, задаются амплитуда (Amplitude), период (Period), ширина (Pulse Width в процентах от периода), фаза (Phase Delay).
Repeating Sequence – последовательность импульсов, их форма задаётся в виде пар чисел (время; величина сигнала).
Sine Wave – синусоидальный сигнал, задаётся амплитуда (Amplitude), частота (Frequency), фаза (Phase) и среднее значение (Bias).
Signal Builder – построитель сигналов, позволяющий задавать форму сигнала, перетаскивая мышью опорные точки.
Random Number – случайные числа с нормальным (гауссовым) распределением. Можно задать среднее значение (Mean Value), дисперсию (Variance), период изменения сигнала (Sample Time).
Uniform Random Number – случайные числа с равномерным распределением в заданном интервале от Minimum до Maximum.
Band Limited White Noise – случайный сигнал, ограниченный по полосе белый шум (имеющий равномерный спектр до некоторой частоты). Блок используется как источник белого шума для моделей непрерывных систем. Задаётся интенсивность (Noise Power) и интервал дискретизации (Sample Time), в течение которого удерживается постоянное значение сигнала. Чем меньше интервал, тем точнее моделирование, однако больше вычислительные затраты.
Display – цифровой дисплей, показывает изменение Scope – осциллограф, показывает изменение сигнала в виде графика, позволяет передавать данные в рабочую область MATLAB для последующей обработки и оформления.
Math Operations Gain – усилитель, задаётся коэффициент усиления (Gain).
Sum – сумматор, используется для сложения и вычитания входов. Параметр List of signs задаёт количество входов, их знаки («+» для сложения и «–» для вычитания). Промежутки между входами (обозначаются знаком |).
Trigonometric Function – тригонометрическая функция.
Signal Routing Manual Switch – ручной переключатель, позволяет двойным щелчком переключать выход на один из двух входных сигналов.
Mux – мультиплексор, объединяет несколько сигналов в один «жгут» (векторный сигнал), в параметрах задаётся Demux – демультиплексор, позволяет «разбить» векторный сигнал на несколько скалярных, в параметрах задаётся число выходов (Number of Outputs).
ИССЛЕДОВАНИЕ ЛИНЕЙНОЙ СТАЦИОНАРНОЙ
ДИНАМИЧЕСКОЙ СИСТЕМЫ В СРЕДЕ MATLAB
Цель работы: освоение методов анализа одномерной линейной непрерывной системы с помощью среды MATLAB.Задачи работы:
• ввести модель системы в виде передаточной функции;
• построить эквивалентные модели в пространстве состояний и в форме «нули-полюса-коэффициент усиления»;
• научиться строить импульсную и переходную характеристики, частотные характеристики;
• научиться использовать окно LTI Viewer для построения различных характеристик;
• научиться строить процессы на выходе линейной системы при произвольном входном сигнале.
Содержание отчёта: исследуемая передаточная функция, копия протокола работы (вводимые команды и результаты вычислений) в среде MATLAB (копируется через буфер обмена из рабочего окна среды MATLAB), графики динамических характеристик.
КРАТКИЕ ТЕОРЕТИЧЕСКИЕ ПОЛОЖЕНИЯ
Модели линейных систем. Для описания линейных стационарных систем могут применяться несколько способов [2]:• дифференциальные уравнения;
• модели в пространстве состояний;
• передаточные функции;
• модели вида «нули-полюса-коэффициенты передачи»;
Первые два способа описывают поведение системы во временной области и отражают внутренние связи между сигналами. Передаточные функции и модели вида «нули-полюса» относятся к частотным способам описания, так как непосредственно связаны с частотными характеристиками системы и отражают только вход-выходные свойства.
Частотные методы позволяют применять для анализа и синтеза алгебраические методы, что часто упрощает расчёты. С другой стороны, для автоматических вычислений более пригодны методы, основанные на моделях в пространстве состояний, поскольку они используют вычислительно устойчивые алгоритмы линейной алгебры.
Исходные уравнения динамики объектов, которые строятся на основе законов физики, в большинстве случаев имеют вид нелинейных дифференциальных уравнений. Для приближённого анализа и синтеза обычно проводят их линеаризацию в окрестности установившегося режима и получают линейные дифференциальные уравнения.
Линеаризованная модель объекта может быть описана дифференциальным уравнением вида:
при нулевых начальных условиях. Причём m n.
Модель в пространстве состояний связана с записью дифференциальных уравнений в стандартной форме Коши (в виде системы уравнений первого порядка):
где x – вектор переменных состояния размера n 1 ; u – вектор входных сигналов (вектор управления) размера m 1 и y – вектор выходных сигналов размера p 1. Кроме того, A, B, C и D – постоянные матрицы.
Согласно правилам матричных вычислений, матрица A должна быть квадратной размера n n, матрица B имеет размер n m, матрица C – p n и матрица D – p m. Для систем с одним входом и одним выходом матрица D – скалярная величина.
Передаточная функция W (s ) линейной стационарной системы от комплексной переменной s определяется как отношение преобразования Лапласа выхода к преобразованию Лапласа входа при нулевых начальных условиях:
Для объекта, описываемого представленным выше дифференциальным уравнением, передаточная функция имеет вид:
По передаточной функции можно построить модель в форме «нулиполюса». Нулями называются корни числителя, полюсами – корни знаменателя.
Динамические характеристики линейных систем. Динамические свойства систем характеризуют реакции на входные воздействия специального вида. Динамические характеристики систем подразделяются на временные и частотные.
Временные характеристики. Импульсной характеристикой (весовой функцией) w(t ) называется реакция системы на единичный бесконечный импульс (дельта-функцию или функцию Дирака) при нулевых начальных условиях. Дельта-функция (t ) определяется равенствами Это обобщённая функция – математический объект, представляющий собой идеальный сигнал, никакое реальное устройство не способно его воспроизвести. Дельта-функцию можно рассматривать как предел прямоугольного импульса единичной площади с центром в точке t = 0 при стремлении ширины импульса к нулю.
Переходной характеристикой (переходной функцией) h(t ) называется реакция системы (при нулевых начальных условиях) на единичный ступенчатый сигнал (единичный скачок):
Импульсная и переходная функции связаны выражениями Для систем без интеграторов переходная характеристика стремится к постоянному значению. Переходная характеристика системы с дифференцирующим звеном (числитель передаточной функции имеет нуль в точке s = 0 ) стремится к нулю. Если система содержит интегрирующие звенья, переходная характеристика асимптотически стремится к прямой, параболе и т.д., в зависимости от количества интеграторов.
По переходной характеристике можно найти важнейшие показатели качества системы – перерегулирование (overshoot) и время переходного процесса (settling time).
Перерегулирование определяется как где hmax – максимальное значение функции h(t ), а h = lim h(t ) – устаt новившееся значение выхода.
Время переходного процесса – это время, после которого сигнал выхода отличается от установившегося значения не более, чем на заданную малую величину (в среде MATLAB по умолчанию используется точность 2%).
Частотные характеристики. Благодаря широкому применению при исследовании устойчивости динамических систем и проектировании регуляторов получили распространение частотные характеристики.
При подаче на вход линейной системы гармонического (синусоидального) сигнала u (t ) = sin ( t ) с частотой (она измеряется в радианах в секунду), на выходе будет также гармонический сигнал той же частоты, но другой амплитуды и фазы y (t ) = A sin( t + ), где A – амплитуда и – сдвиг фазы.
Для построения частотной характеристики надо использовать подстановку s = j в передаточной функции W (s). Выражение W ( j) называется частотной передаточной функцией или амплитудно-фазовой частотной характеристикой системы (АФЧХ).
Зависимость модуля величины W ( j) от частоты называется амплитудной частотной характеристикой (АЧХ), а зависимость аргумента комплексного числа (фазы) W ( j) от частоты – фазовой частотной характеристикой (ФЧХ):
АЧХ показывает, насколько усиливается амплитуда сигналов разных частот после прохождения через систему, а ФЧХ характеризует сдвиг фазы сигнала.
На практике широкое применение получила диаграмма Боде (логарифмическая амплитудная характеристика – ЛАХ), которая определяется как L() = 20 lg A( w). Измеряется ЛАХ в децибелах и строится как функция от lg(w).
Пакет CONTROL System Toolbox. В состав системы MATLAB входит пакет прикладных программ (ППП) CONTROL System Toolbox, который предназначен для работы LTI-моделями (Linear Time Invariant Models – линейные инвариантные во времени системы) систем управления.
В русскоязычной литературе за такими системами закрепилось название – линейные стационарные системы.
В состав этого пакета входит более ста различных функций, в том числе функции создания и преобразования моделей, представленных в виде уравнений состояния, передаточных функций, функций для построения частотных и временных характеристик системы и др.
Для построения передаточной функции используется функция tf.
Она имеет вид:
где В – вектор коэффициентов числителя передаточной функции, записанных по убыванию степеней; А – вектор коэффициентов знаменателя передаточной функции.
П р и м е р : необходимо сформировать передаточную функцию В данном случае векторы коэффициентов числителя и знаменателя передаточной функции имеют вид: B=[2 4], A=[1 0.8 2.06 2.26].
Процедура образования передаточной функции в диалоговом режиме имеет вид:
» B=[2 4];
» A=[1 0.8 2.06 2.26];
» f=tf(B,A) После нажатия клавиши в памяти создаётся объект класса tf, описывающий передаточную функцию, и на экране появится передаточная функция в виде:
Transfer function:
----------------------------s^3 + 0.8 s^2 + 2.06 s + 2. Определение передаточной функции F (s ) можно также осуществить без предварительного построения числителя и знаменателя:
» f=tf([2 4], [1 0.8 2.06 2.26]) По передаточной функции можно легко построить модель в форме «нули-полюса-коэффициент передачи»:
» f_zpk = zpk(f) Zero/pole/gain:
-----------------------------s+1) (s^2 – 0.2s + 2.26) Нулями называются корни числителя, полюсами – корни знаменателя. Эта функция имеет один нуль в точке s = 2 и три полюса в точках s = 1 и s = 0,1 ± 1,5i. Паре комплексных полюсов соответствует квадратный трёхчлен s 2 0,2s + 2,26.
Для определения нулей и полюсов передаточной функции F (s ) используются соответственно функции zero и pole:
» zero(f) ans = » pole(f) 0.1000 + 1.5000i 0.1000 – 1.5000i –1. Для преобразования передаточной функции F (s ) в модель в пространстве состояний используется команда » f_ss = ss(f) Continuous-time model.
Это означает, что матрицы модели имеют вид На практике достаточно часто реальные системы автоматического управления (САУ) представляются структурно в виде соединённых между собой отдельных блоков (динамических звеньев), уравнения которых обычно достаточно просты.
В MATLAB предусмотрена возможность программно строить схему САУ путём предварительного ввода моделей звеньев, составляющих САУ, и последующего соединения этих звеньев в единую структуру.
К процедурам, осуществляющим расчёт характеристик соединений отдельных звеньев, относятся:
• plus (minus) или ' + ' (' – ') – выполняет параллельное соединение указанных при обращении звеньев, т.е. определяет характеристики модели системы из параллельно соединённых звеньев, например W1 и W2:
• parallel – осуществляет ту же процедуру параллельного соединения звеньев; в отличие от предыдущей процедуры может использоваться для многомерных систем и осуществления параллельного соединения лишь по некоторым входам и выходам;
• mtimes или ‘ * ’ – осуществляет последовательное соединение звеньев, имена которых указаны (применяется для одномерных систем):
• series – последовательное соединение (может использоваться для многомерных систем);
• feedback – соединение с обратной связью, т.е. такое соединение двух звеньев, когда второе указанное звено составляет цепь отрицательной обратной связи для первого звена:
Пример: Разработать модель замкнутой системы автоматического управления, в которой объект управления представляется в виде соединения трёх апериодических звеньев 1 порядка с передаточными функциями Решение:
»w1=tf([2],[3 1]) »w2=tf([1],[2 1]) »w3=tf([1],[4 1]) »w=w1*w2*w »R=tf([0.5],[1])+tf([0.2],[1 0]) »W_zs= feedback(w,R) »step(W_zs) Пакет CONTROL System Toolbox предоставляет возможности по построению временных и частотных характеристик систем с использованием следующих функций (табл. 1):
step() nyquist() Другим вариантом получения графиков динамических характеристик системы является использование графического окна – LTI Viewer. Вызов его осуществляется командой ltiview, которой, в качестве параметра, можно указать имя переменной, содержащей LTI-объект, например:
» ltiview(f) Загрузку исследуемого объекта можно осуществить также непосредственно из окна LTI Viewer, вызвав меню File, а затем команду Import… Используя команду Plot Configuration…, из меню Edit можно выбрать конфигурацию области отображения графиков (количество и расположение окон построения графиков), а также виды выстраиваемых характеристик системы (Step, Impulse, Bode и др.) – раздел Response type (рис. 1).
Чтобы построить частотную характеристику в MATLAB можно использовать функцию freqresp(f,w), где f – модель линейной системы, заданная в виде передаточной функции, в пространстве состояний или в форме «нули-полюса-коэффициент усиления», w – вектор частот.
Предварительно до использования этой функции необходимо создать массив частот в нужном диапазоне. Для этой цели могут быть задействованы функции linspace и logspace. Первая функция осуществляет равномерное распределение точек по линейной шкале, вторая – по логарифмической.
Процедура построения АЧХ системы на участке частот [0, 10] в диалоговом режиме имеет вид:
» w= linspace (0, 10, 100);
» freq = freqresp(f, w);
» freq = freq(:);
» plot (w, abs(freq));
Для вычисления фазы (в градусах) и построения её на графике могут использоваться команды:
» phi = angle(freq)*180/pi;
» plot (w, phi);
ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ
1. Запустить среду MATLAB и в командном режиме очистить рабочее пространство и окно:» clear all » clc 2. Введите передаточную функцию как объект tf (коэффициенты полиномов передаточной функции вводятся в соответствии с вариантом (табл. 2), определяемым преподавателем):
» A = [a4 a3 a2 a1 a0] » B = [b2 b1 b0] » f = tf(B,A) 3. Найдите нули, полюса передаточной функции и коэффициент усиления:
» z = zero(f) » p = pole(f) 4. Постройте модель исходной системы в форме «нули-полюса»:
» f_zpk = zpk(f) » k = dcgain(f) 5. Постройте модель системы в пространстве состояния:
» f_ss = ss(f) 6. В командном режиме постройте переходную характеристику:
» step(f) 7. В командном режиме постройте импульсную характеристику:
» impulse(f) 8. Постройте логарифмические частотные характеристики (диаграммы Боде):
» bode(f) 9. Постройте амплитудно-фазовую частотную характеристику (частотного годографа Найквиста):
» nyquist(f) 10. Запустите графическое окно LTI Viewer:
» ltiview 11. Загрузите объект f. Для этого в меню File необходимо выбрать пункт Import…, а далее выбрать объект f.
12. Постройте вышеперечисленные динамические характеристики, используя интерфейс LTI Viewer. По окончании работы с LTI Viewer закройте все окна за исключением командного окна MATLAB.
13. Постройте сигнал, имитирующий прямоугольные импульсы единичной амплитуды с периодом 10 секунд:
» [u,t] = gensig('square',10);
14. Выполните моделирование и постройте на графике сигнал выхода системы f при данном входном сигнале u(t):
» lsim (f, u, t)
КОНТРОЛЬНЫЕ ВОПРОСЫ
1. Какие Вы знаете способы описания линейных стационарных динамических систем?2. Задано дифференциальное уравнение системы:
Представьте описание системы в виде передаточной функции; модели вида «нули-полюса»; уравнений состояния.
3. Представьте исследуемую вами систему в виде последовательного соединения типовых звеньев.
4. Постройте динамические характеристики типовых звеньев.
5. Постройте динамические характеристики типовых регуляторов (П, И, ПИ, ПИД, ПД).
№ Вид передаточной функции
МОДЕЛИРОВАНИЕ СИСТЕМ УПРАВЛЕНИЯ В
ПАКЕТЕ SIMULINK
Цель работы: освоение методов моделирования линейных систем в пакете SIMULINK Задачи работы:• научиться строить и редактировать модели систем управления в пакете SIMULINK;
• научиться изменять параметры блоков;
• научиться строить переходные процессы;
• оценить влияние настроечных параметров ПИД-регулятора на качественные показатели процесса регулирования в одноконтурной АСР.
Содержание отчёта: краткое описание исследуемой системы, схемы ПИД-регулятора и одноконтурной АСР из пакета SIMULINK, графики переходных процессов в одноконтурной АСР при изменении задания и отработке возмущения с различными типами регуляторов и соответствующие им настройки регуляторов, анализ влияния настроечных параметров на процесс регулирования.
ОПИСАНИЕ МОДЕЛИРУЕМОЙ СИСТЕМЫ
В работе требуется провести исследование одноконтурной АСР с ПИД-регулятором. Её структурная схема показана на рис. 2.Рис. 2. Структурная схема одноконтурной АСР Передаточная функция промышленных объектов во многих случаях с достаточной точностью может быть представлена в виде:
где T, – соответственно, большая и меньшая постоянные времени объекта управления; K об – коэффициент усиления объекта управления; – время чистого запаздывания.
Передаточная функция ПИД-регулятора где K p, K и, K д – соответственно, настроечные коэффициенты пропорциональной, интегральной и дифференциальной составляющей.
ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ
1. Для запуска пакета SIMULINK щёлкните по кнопке в командном окне MATLAB или введите команду simulink в командной строке.2. Создайте новую модель с помощью верхнего меню открывшегося окна Simulink Library Browser. Для этого в меню File выберите пункт New, а затем Model, или воспользуйтесь сочетанием клавиш Cntl-N.
3. Перетащите блок Gain (усилитель) из окна Simulink Library Browser (группа Math Operations) в окно модели и дважды продублируйте его. Отредактируйте описание блоков – Kp, Ki, Kd.
4. Перетащите блоки Integrator (интегратор) и дифференциатор из окна Simulink Library Browser (группа Continuous) в окно модели.
5. Перетащите блок Sum (сумматор) из Simulink Library Browser (группа Math Operations) в окно модели и перейдите в окно редактирования его свойств (дважды нажав на левую клавишу мыши на изображении блока). Для свойства Icon shape установите rectangular (прямоугольный вид сумматора), а для свойства List of signs задайте +++ (три плюса), говорящие о том, что будет использовано три суммирующих входа.
6. Из Simulink Library Browser, соответственно из групп Sources и Sinks перетащить блоки In1 и Out1.
7. Произведите соединение блоков в соответствии со схемой, представленной на рис. 3. Схема реализует ПИД-регулятор (1).
8. Создайте из нарисованной схемы, реализующей ПИД-регулятор, подсистему (Subsystem). Для этого выделите все элементы схемы, нажмите на выделенной области правую клавишу мыши и в появившемся окне выберите пункт Create Subsystem.
9. Отредактируйте название образовавшегося блока, назвав его PID.
10. Перетащите 2 блока Transfer Fcn (передаточная фунция) и Transport Delay (Запаздывание) из окна Simulink Library Browser (группа Continuous) в окно модели. Введите числитель и знаменатель передаточной функции в соответствии с номером варианта (табл. 3), а также время запаздывания Time Delay для блока Transport Delay.
Denominator [T 1] 11. Перетащите 2 блока Sum (сумматор) из Simulink Library Browser в окно модели. Для организации отрицательной обратной связи отредактируйте свойство List of signs (задайте + –) одного из блоков.
12. Перетащите в окно модели блок Scope (осциллограф) из группы Sinks и установите его в правой части окна.
13. Перетащите в окно модели 2 блока Step из группы Sources и установите его слева от сумматора. Дайте название первому блоку Zadanie, второму – Vozmuschenie.
14. Установите в обоих блоках время подачи сигнала равное 0.
15. Произведите соединение блоков в соответствии со схемой, представленной на рис. 4. Схема реализует одноконтурную АСР.
16. Сохраните модель в своей папке под именем lab2.mdl 17. Установите время моделирования 100 секунд.
18. В блоке с именем Zadanie (тип блока Step) установите величину подаваемого сигнала равную 1, а в блоке Vozmuschenie – 0. При этом при выполнении моделирования будет имитироваться ситуация отработки изменения задания на 1. При отработке возмущения необходимо наоборот установить 0 для блока Zadanie и 1 – для блока Vozmuschenie.
Vozmuschenie Рис. 4. Реализация одноконтурной АСР в пакете SIMULINK 19. Выполните моделирование, нажав левой клавишей мыши на кнопке. Вид переходного процесса в одноконтурной АСР при изменении задания можно наблюдать в окне Scope.
20. Подберите настроечные параметры регулятора, обеспечивающие устойчивость замкнутой системы регулирования. Настройки регуляторов (Kp, Ki, Kd) устанавливаются при редактировании блоков Gain в блоке PID.
21. Постройте график переходного процесса при отработке возмущения.
22. Проведите исследование одноконтурной АСР с различными типами регуляторов (П, ПИ, ПИД, ПД) с целью оценки влияния настроечных параметров регулятора на прямые показатели качества процесса регулирования (статическая и динамическая ошибки регулирования, величина перерегулирования, время переходного процесса), а также на устойчивость системы.
КОНТРОЛЬНЫЕ ВОПРОСЫ
1. Как запустить пакет SIMULINK?2. Что такое Simulink Library Browser?
3. Какое расширение имеют файлы – модели SIMULINK?
4. Как создать новую модель?
5. Как создать подсистему (Subsystem)?
6. Как сделать, чтобы один и тот же сигнал поступал на несколько блоков?
7. Как скопировать (продублировать) блок в окне модели?
8. Как скопировать изображение модели в документ Microsoft Word?
9. Как изменить знаки арифметических действий в сумматоре?
10. Как ввести параметры блока Transfer Fcn (передаточная функция)?
11. Как изменить время моделирования?
12. Почему при использовании П- или ПД-регулятора система не компенсирует постоянное возмущение?
13. Какими свойствами должен обладать регулятор для того, чтобы постоянное возмущение полностью компенсировалось?
14. С помощью каких показателей качества оценивают переходный процесс в АСР?
ПОСТРОЕНИЕ МАТЕМАТИЧЕСКОГО ОПИСАНИЯ ОБЪЕКТА
УПРАВЛЕНИЯ ЭКСПЕРИМЕНТАЛЬНЫМИ МЕТОДАМИ
Цели работы: освоение методики построения математического описания объекта регулирования по результатам экспериментально полученных статических и динамических характеристик объекта.Задачи работы:
• освоение методики аппроксимации экспериментальной статической характеристики объекта регулирования;
• освоение одной из методик аппроксимации экспериментальной динамической характеристики объекта регулирования передаточной функцией;
• закрепление навыков по созданию и использованию m-файлов;
• знакомство с методикой решения оптимизационных задач.
Содержание отчёта: исходные данные, результаты аппроксимации статической характеристики (аппроксимирующий полином, графики), результаты аппроксимации динамической характеристики (передаточная функция, графики).
КРАТКИЕ ТЕОРЕТИЧЕСКИЕ ПОЛОЖЕНИЯ
Математические модели объектов управления, процедура их построения зависят от их целевого назначения, свойств объекта, а также функций и задач системы управления.При автоматизации производственных процессов наиболее правильные решения могут быть приняты на основании данных о статических и динамических свойствах регулируемого объекта. Для получения этих данных определяют соответственно статические и динамические характеристики регулируемых объектов.
Статические характеристики АСР или отдельных её звеньев представляют собой функциональную зависимость выходных величин от входных при статическом (стационарном, установившемся) режиме работы, т.е. при не изменяющихся во времени значениях входных и выходных величин.
Статические характеристики позволяют определить значения регулируемых параметров, положение регулирующих органов, расходы веществ или энергии и другие данные для любого установившегося состояния исследуемого объекта. Эти характеристики используют при расчётах АСР, расчётах оптимальных режимов, при оценке линейности объектов, при технологических и экономических расчётах.
Динамические характеристики представляют собой зависимости между изменениями входных и выходных величин в динамическом режиме (во времени). Они могут быть представлены в виде дифференциальных уравнений, передаточных функций, временных и частотных характеристик. Имея динамическую характеристику, представленную в каком-либо виде, нетрудно преобразовать её и представить в любом другом (из перечисленных выше) виде.
Динамические характеристики дают информацию об инерционных свойствах регулируемых объектов (систем, элементов систем) и таким образом являются исходными данными при синтезе автоматических систем регулирования. Они позволяют выполнить эту работу в полном объёме и завершить её расчётом параметров настройки на данном регулируемом объекте с целью получения желаемых форм графиков процессов регулирования.
Статические и динамические характеристики могут быть получены аналитическим, экспериментальным или экспериментально-аналитическим методом.
При построении систем автоматической стабилизации отдельных технологических параметров достаточно располагать динамическими и статическими характеристиками объекта управления в узком диапазоне изменения входных и выходных координат. Для получения этих характеристик применяются многочисленные экспериментальные методы. Эти методы базируются на предположении о линейности и сосредоточенности параметров объекта, неизменности во времени его статических и динамических характеристик. Принятие этих допущений позволяет достаточно просто описать изменение наблюдаемых выходных координат объекта линейными дифференциальными уравнениями с постоянными коэффициентами.
При построении экспериментальных математических моделей используют пассивный или активный эксперименты.
При активных методах на входе исследуемого объекта создаются искусственные испытательные воздействия (ступенчатые, импульсные, гармонические и др.). Пассивными методами называют такие, при которых не допускаются искусственные возмущения на входе исследуемого объекта, а в качестве входных испытательных воздействий используются естественные воздействия, возникающие в процессе нормальной эксплуатации.
Аппроксимация статических характеристик. При экспериментальном исследовании статических характеристик некоторого одноканального объекта управления получено m значений входной xi i = 1, m и выходной yi координат. Требуется найти функциональную зависимость y = f ( x, a), наилучшим образом описывающую экспериментальные данные, где a – вектор параметров функции. Для этой цели достаточно часто используется метод наименьших квадратов. Согласно этому методу вводится функция, характеризующая степень близости экспериментальных и расчётных данных:
При использовании метода наименьших квадратов задаются видом функции и вектором параметров этой функции. Часто в качестве аппроксимирующей функции y = f ( x, a) используют степенную функцию вида:
Задача поиска аппроксимирующей функции (вектора параметров) ставится и решается как задача безусловной многомерной оптимизации: требуется найти вектор параметров a, обеспечивающий минимум функции S.
При использовании пакета MATLAB полиномиальная аппроксимация данных измерений, которые сформированы как некоторый вектор Y, при некоторых значениях аргумента, которые образуют вектор Х такой же длины, что и вектор Y, осуществляется процедурой polyfit(X, Y, n). Здесь n – порядок аппроксимирующего полинома. Результатом действия этой процедуры является вектор длиной (n + 1) из коэффициентов аппроксимирующего полинома.
Аппроксимация динамических характеристик. В инженерной практике при настройке АСР математическое описание объекта регулирования достаточно часто представляется в виде передаточных функций.
Передаточная функция объекта регулирования может быть определена в результате аппроксимации экспериментальной динамической характеристики, полученной в результате проведения активного эксперимента.
В качестве испытательного воздействия чаще всего используется ступенчатое воздействие.
Алгоритм определения аппроксимирующей передаточной функции объекта регулирования имеет вид:
• При подаче на вход объекта ступенчатого воздействия на выходе снимается переходная характеристика h(t ).
• По виду кривой определяется структура передаточной функции.
Для большого числа реальных объектов регулирования модель объекта может быть представлена в виде:
• Аппроксимация переходной характеристики объекта передаточной функцией. Аппроксимация заключается в подборе таких параметров передаточной функции, которые обеспечивают максимальную близость экспериментальной кривой разгона h(t ) и кривой разгона, соответствующей передаточной функции (2) – hp (t ).
Настроечные параметры k, T1, T2, передаточной функции (2) могут быть получены при минимизации функции
ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ
Аппроксимация статической характеристики 1. Запустите среду MATLAB. В командном режиме очистите рабочее пространство и окно:» clear all » clc 2. Сформируйте массивы экспериментальных точек – векторы Х и Y (в соответствии с вариантом (табл.4), определяемым преподавателем).
» Х = [х1 х2 х3 х4 х5 х6 х7 х8 х9 x10];
» Y = [y1 y2 y3 y4 y5 y6 y7 y8 y9 y10];
3. Произведите аппроксимацию экспериментальных точек полиномом первой степени (линейной функцией) и второй степени и определите векторы коэффициентов полиномов – соответственно P1 и P2.
» P1 = polyfit(X,Y,1) » P2 = polyfit(X,Y,2) 4. Используя функцию polyval, постройте аппроксимирующие полиномы.
» F1 = polyval(P1,X) » F2 = polyval(P2,X) 5. На одном рисунке постройте экспериментальные точки, а также аппроксимирующие полиномы первой и второй степени » plot(X,Y,'o',X,F1,'-.',X,F2,'-') Аппроксимация динамической характеристики Исходными данными для аппроксимации являются данные об экспериментальной переходной характеристике, представленной в табличной форме, например:
h 0 0 0 0.07 0.23 0.43 0.64 0.84 1.03 1.19 1.33 1.45 1.55 1.64 1.71 1. Вариант исходных данных определяется преподавателем (табл. 5).
4. Создайте Script-файл для решения задачи (3). Для этого вызовите меню File командного окна MATLAB и выберите в нём сначала команду New, а затем команду M-file.
5. В появившемся окне текстового редактора наберите текст:
% t – Массив моментов времени % h – массив значений экспериментальной переходной характеристики % hp – массив значений рассчитанной переходной характеристики % a – вектор параметров передаточной функции объекта % a(1)=k, a(2)=t1, a(3)=t2, a(4)=tau global t h hp t=0:1:15;
h=[0 0 0 0.07 0.23 0.43 0.64 0.84 1.03 1.19 1.33 1.45 1.55 1.64 1.71 1.77];
[a,fval] = fminunc(@func,[1 1 2 1]) plot(t,h,'o',t,hp,'-') Процедура fminunc осуществляет поиск минимума функции нескольких переменных. Имя минимизируемой функции задаётся в первом аргументе процедуры fminunc, а во втором – передаётся начальная точка для поиска.
6. Сохраните этот текст в файле под именем din_app.m.
7. Создайте файл-функцию, высчитывающую минимизируемую функцию (3). Текст файл-функции имеет вид:
function f=func(a) global t h hp y=series(tf([a(1)],[a(2) 1]),tf([1],[a(3) 1]));
y.inputd = a(4);
[hp,t]=step(y,t);
f=sum((h-hp').^2);
В этой функции строки y=series(tf([a(1)],[a(2) 1]),tf([1],[a(3) 1]));
y.inputd = a(4);
определяют передаточную функцию объекта в соответствии с (2). Здесь a(1), a(2), a(3), a(4) соответственно k, T1, T2,.
Строка [hp,t]=step(y,t); строит отклик на ступенчатое воздействие.
Последняя строка определяет сумму квадратов невязок между экспериментальной и рассчитанной переходными характеристиками.
8. Сохраните введённый текст в файле под именем func.m.
9. В командном режиме запустите на исполнение Script-файл:
» din_app 10. При завершении выполнения этого Script-файла в командном окне MATLAB будет отображён результат решения задачи минимизации функции (3) в виде:
2.0382 2.4265 4.8308 2. fval = 3.9104e– Также на графике (рис. 5) будут представлены экспериментально снятые точки переходной характеристики и рассчитанная аппроксимирующая переходная характеристика.
Рис. 5. Результаты аппроксимации переходной характеристики 4. Таблица вариантов экспериментальных статических 5. Таблица вариантов экспериментальных динамических характеристик
КОНТРОЛЬНЫЕ ВОПРОСЫ
1. Что позволяют определить статические характеристики объекта?2. Как осуществить линеаризацию статической характеристики?
3. Какие виды динамических характеристик вы знаете?
4. Какая взаимосвязь между импульсной и переходной характеристиками?
5. Какие методы аппроксимации динамических характеристик вы знаете?
6. Чем отличается файл-функция от Script-файла в пакете MATLAB?
7. Какая функция может быть использована для проведения аппроксимации полиномом n-й степени?
8. Что осуществляет функция polyval в пакете MATLAB?
9. Какая функция в пакете MATLAB может быть использована для поиска минимума функции нескольких переменных?
ПОСТРОЕНИЕ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ
АНАЛИТИЧЕСКИМ МЕТОДОМ
(МОДЕЛИРОВАНИЕ ТЕПЛООБМЕННОЙ АППАРАТУРЫ)
Цель работы: ознакомление с методами моделирования процесса теплообмена через стенку и с расчётом теплообменных аппаратов на ЭВМ.Задания работы:
1. Ознакомиться с принципами составления тепловых балансов потоков и разделяющей стенки в теплообменных и реакционных аппаратах химической технологии при теплообмене через стенку.
2. Отметить особенности математических описаний при разных способах организации теплообмена.
3. Выполнить вариант задания на ЭВМ.
Содержание отчёта:
1. Вариант задания.
2. Математическая модель теплообменного аппарата.
3. Все расчётные формулы, используемые в работе.
4. Графики рассчитанных зависимостей для теплоносителей.
5. Тексты программ.
КРАТКИЕ ТЕОРЕТИЧЕСКИЕ ПОЛОЖЕНИЯ
Процесс передачи тепла через стенку весьма распространён в химической технологии и значительно влияет на протекание химических реакции во всех типах реакторов. Процесс передачи тепла в теплообменной аппаратуре является основным и служит для сообщения технологическому потоку нужной температуры.Выбирая различные способы оформления реакторов, можно влиять на интенсивность теплообмена между основным (реакционным) потоком и потоком хладагента или окружающей средой. При полном отсутствии теплообмена через стенку получают адиабатический реактор. Реакторы, имеющие теплообмен с внешней средой, относятся к политропическим.
При рассмотрении процесса передачи тепла от одного теплоносителя к другому через стенку можно выделить несколько элементарных этапов:
переход тепла от горячего теплоносителя к более холодной стенке, поглощение тепла материалом стенки и её нагрев, распределение тепла по объёму стенки, переход тепла от стенки к холодному теплоносителю.
Если процесс теплообмена протекает стационарно, то температура в каждой точке материала (теплоносителей и стенки) не изменяется во времени. Применение модели с сосредоточенными параметрами (т.е. когда пространственные координаты не входят в математическое описание) приводит к алгебраическим соотношениям между температурами в системе. Если, наоборот, температуры меняются во времени, математическое описание получается в виде системы обыкновенных дифференциальных уравнений (аргументом является время).
Зависимость температур от геометрических координат обуславливает математическое описание статики в виде обыкновенных дифференциальных уравнений (если пространственная координата одна) или дифференциальных уравнений в частных производных. Независимыми переменными при этом являются пространственные координаты. Динамическая модель при наличии пространственно-распределённых эффектов описывается уравнениями в частных производных, причём одной из независимых переменных является время.
Интенсивность перехода тепла от одного теплоносителя (например, горячего потока жидкости или газа) к другому (стенке) зависит от разности температур между ними, а также от теплового сопротивления. В расчётные уравнения, однако, обычно включают не сопротивление, а обратную величину – коэффициент теплоотдачи, q – тепловой поток (ккал/ч или Вт) через поверхность площадью 1 м2 при разности температур (температурном напоре) 1 °С.
Полный тепловой поток q определяется произведением коэффициента теплоотдачи, площади поверхности F и температурного напора T:
Уравнение (4) применимо как к нагреванию стенки от горячей жидкости, так и, наоборот, к нагреванию холодной жидкости горячей стенкой;
при этом T будет иметь разные знаки.
Если пренебречь распространением тепла в стенке, то теплопередачу от горячего потока жидкости к холодному, находящемуся по другую сторону стенки, можно представить как процесс преодоления тепловым потоком двух последовательных сопротивлений теплоотдачи – от горячего потока к стенке и от нагревшейся стенки к холодному потоку.
Используя вместо сопротивлений коэффициенты теплоотдачи (1 и 2), получаем выражение, определяющее коэффициент теплопередачи (K):
В практических расчётах часто используют коэффициент теплопередачи как характеристику интенсивности теплообмена между потоками:
В тех случаях, когда коэффициенты теплоотдачи учитываются порознь, принимают усреднённую температуру стенки, разделяющей потоки. Иными словами, считают, что теплопроводность материала стенки настолько велика, что перепад температуры отсутствует.
Коэффициенты теплоотдачи зависят от многих параметров, но наиболее сильно – от скорости потока, характера набегания жидкости на стенку, плотности и теплопроводности жидкости. При выполнении точных расчётов зависимость коэффициента теплоотдачи от параметров потока следует учитывать. Однако для большинства инженерных расчётов теплообменной аппаратуры и реакторов достаточны упрощённые представления.
Для вывода уравнений математического описания процесса теплообмена через стенку следует рассмотреть тепловой баланс каждой среды, имеющей запас тепла. Он складывается из прихода и расхода тепла, которые определяют накопление тепла в объёме; накопление является временным процессом: накопление = приход – расход. В статике ввиду равенства прихода и расхода тепла накопление равно нулю.
Накопление связано с изменением температуры:
или для элементарного объёма где – плотность; сp – удельная теплоёмкость; V – объём; S – сечение потока; d1 – элементарный участок потока.
Приход и расход тепла могут определяться теплоотдачей (теплопередачей), а в случае проточной системы с распределёнными параметрами – притоком и уносом тепла с конвективным потоком.
Количество тепла, поступающее в аппарат с конвективным потоком, определяется как или для элементарного объёма за элементарное время d где G – объёмный расход потока.
Количество тепла, уходящее из рассматриваемого объёма с конвективным потоком, определяется следующим выражением:
или для элементарного объёма Приход тепла, определяемый теплопередачей:
или для элементарного объёма за элементарное время где Твн – температура внешнего теплоносителя.
С учётом полученных соотношений накопление тепла в системе составит:
или в элементарном объёме за элементарное время Проведя несложные преобразования, получим уравнение теплового баланса, описывающее динамику теплообменников, во всём объёме которых происходит полное (идеальное) смешение частиц потока:
где Т0, Т – температура потока на входе и в зоне идеального смешения.
Соответственно для трубчатых теплообменников, работающих по принципу вытеснения, уравнение динамики будет выглядеть следующим образом:
Ввиду того, что в статическом режиме накопление тепла в системе равно нулю, модель статики теплообменников смешения будет иметь вид:
статика трубчатых теплообменников описывается уравнением где u – скорость движения потока.
Пример 1. Теплообменник представляет собой тонкостенный змеевик, по которому в режиме идеального вытеснения движется охлаждаемый поток жидкости. Змеевик погружён в воду, непрерывно протекающую через сосуд, так что температура охлаждающей воды Твн практически постоянна и равна 10 °С во всём объёме.
Требуется определить температуру на выходе потока, идущего по змеевику со скоростью u = 4 м/с, если температура его на входе равна 95 °С, длина трубки змеевика L = 2 м, его сечение S = 10–4 м2, коэффициент теплопередачи K = 1,16·104 Вт/(м2 °С). Теплоёмкость охлаждаемой жидкости cp = 2,93·103 Дж/(кг °С), её плотность = 900 кг/м3. Параметры считать не зависящими от температуры; изменение объёма не учитывать.
Режим работы считать стационарным.
Температура охлаждаемого потока Т подчиняется дифференциальному уравнению (10):
где l – длина; r – радиус змеевика; uS – объёмный расход потока.
Начальное условие для уравнения (11) Т(0) = 95 °С.
Вычислим коэффициент уравнения:
Уравнение (12) подлежит решению в пределах изменения независимой переменной 1 от 0 до 2 м.
Решение представлено на рис. 6. Температура на выходе потока составила T(L) = 47 °С.
Для численного интегрирования уравнения (12) в MATLAB необходимо создать два m-файла. В файле func1_T.m описывается функция для вычисления правой части дифференциального уравнения (12). В файле Tepl1.m осуществляются задание исходных данных для решения задачи и начального условия, а также решение дифференциального уравнения с использованием функции Matlab ode45.
Рис. 6. Профиль температуры по длине теплообменника Функция ode45 в данном случае имеет следующий синтаксис:
где @func1_T – дескриптор ODE-функции правой части дифференциального уравнения; [l, T] – матрица решений T, где каждая строка соответствует длине, возвращённой в векторе-столбце l; [0 L] – вектор, определяющий интервал интегрирования; T0 – вектор начальных условий.
function dT = func1_T(l,T) % Функция правой части Tvn=10;
dT = 0.39*(Tvn-T);
L=2;
T0=95;
[l,T] = ode45(@func1_T,[0 L],T0); % Функция решения plot(l,T);
Пример 2. Жидкость охлаждается в теплообменнике типа «труба в трубе». Охлаждаемая жидкость и хладагент движутся параллельно (прямотоком). Требуется определить температуры потоков на выходе теплообменника, если начальная температура охлаждаемой жидкости равна 170 °С, а хладагента 15 °С. Плотность охлаждаемой жидкости и хладагента = 900 кг/м3. Диаметры труб теплообменника: внутренней D1 = 0,1 м, наружной (для хладагента) D2 = 0,3 м. Длина теплообменника L = 1 м. Тепломкость жидкости и хладагента сp = 3,35·103 Дж/(кг °С). Объёмный расход охлаждаемой жидкости G1 = 2,28·10–4 м3/с, хладагента G2 = 5,75·10-4 м3/с, коэффициент теплопередачи K = 4900 Вт/(м2°С). Температурный профиль по длине для каждого из потоков определяется решением системы дифференциальных уравнений:
где T1 и Т2 – температуры охлаждаемой и охлаждающей жидкости.
Начальные условия: Т1(0) = 170 °С; Т2(0) = 15 °С.
После подстановки в уравнения (13) и (14) численных значений параметров получаем следующую систему:
Графики решения системы уравнений математического описания статики теплообменника представлены на рис. 7. На нём изображены температурные профили вдоль теплообменника для обоих теплоносителей.
Можно видеть, что движущая сила процесса сильно меняется по длине, поэтому эффективность использования различных участков теплообменника не одинакова. Температуры теплоносителей на выходе теплообменника равны: T1(L) = 62,9 °С, Т2(L) = 56,5 °С.
Численное интегрирование уравнений (15) и (16) аналогично примеру 1. Создаются два m-файла: func2_T.m – функция для вычисления правых частей дифференциальных уравнений, Tepl1.m – задание исходных данных, начальных условий и решение дифференциальных уравнений с использованием функции MATLAB ode45.
Рис. 7. Изменение температур теплоносителей Файл func2_T.m function dT = func2_T(l,T) % Функция правых частей dT = zeros(2,1);
dT(1)=2.239*(T(2)-T(1)); % Уравнение правой части 1-го dT(2)=0.888*(T(1)-T(2)); % Уравнение правой части 2-го 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'); % Построение зависимостей Пример 3. Смоделировать статический режим теплообменника типа «труба в трубе», используя данные, приведённые в примере 2, для случая противотока. Принять полную длину теплообменника L = 2,5 м.
Тепловые процессы в противоточном теплообменнике подчиняются тем же закономерностям, что и в прямоточном. Поэтому математическое описание теплообменника записывается аналогично, однако формально однотипные уравнения для обоих теплоносителей имеют аргументы различного знака:
Существенным различием, отражающим иную организацию потоков теплоносителей, является принципиально другое задание условий решения уравнений (17) и (18) по сравнению с заданием при решении уравнений (13) и (14).
Совместное интегрирование уравнений (17) и (18) возможно лишь в одном направлении: либо при l, меняющемся от 0 до L, либо в обратном – от L до 0. При этом в любом случае оговорено лишь одно начальное условие, второе остаётся неизвестным. Известно лишь, к какому значению в конце решения должна подойти вторая переменная.
Для решения задачи воспользуемся последним обстоятельством и попытаемся отыскать неизвестное начальное условие Т2(0) с таким расчётом, чтобы условие, заданное для конца решения (граничное условие), было выполнено, т.е. T2(L) = 15 °С. Такие задачи при малом числе условий, подлежащих определению, обычно решают методом проб и ошибок.
Задачей поиска начального условия Т2(0) является выполнение граничного условия T2(L) при интегрировании системы уравнений (17) и (18).
Рисунок 8, а иллюстрирует процесс поиска неизвестного начального условия Т2(0). Кривая 1 отражает профиль температуры Т2, полученный в предположении, что хладагент нагреется до 57 °С (значение взято произвольно). Видно, что это предположение неудачно, так как не получено ожидаемое значение T2(L) = 15 °С.
Рис. 8. Результаты моделирования на ЭВМ противоточного теплообменника:
а – серия кривых решения системы уравнений (17) и (18), полученная в процессе определения недостающего начального условия;
б – изменение температур теплоносителей по длине теплообменника Приняв Т2(0) = 65 °С, получаем кривую 2. Она также неудовлетворительна, так как вместо Т2(L) = 15 °С получаем лишь 4 °С. Кривые 3 и также не соответствуют искомому начальному условию: в первом случае входная температура T2(L) получилась заниженной, во втором – завышенной. По-видимому, искомое начальное условие лежит между двумя последними проверенными значениями (70 и 80 °С).
Тщательно исследуя намеченный диапазон задания начального условия, находим Т2(0) = 75 °С.
На рисунке 8, б показан результат решения задачи. Выходная температура охлаждаемого потока T1(L) достигает 18 °С.
Численное интегрирование уравнений (17) и (18) аналогично примеру 2.
Файл func3_T.m function dT = func3_T(l,T) % Функция правых частей global b1 b2;
dT = zeros(2,1);
dT(1)=b1*(T(2)-T(1)); % Уравнение правой части 1-го dT(2)=-b2*(T(1)-T(2)); % Уравнение правой части 2-го 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');
Пример 4. Смоделировать переходный режим теплообменника типа «смешение – смешение». Теплообменник представляет собой двухкамерную ёмкость. В первую камеру ёмкости поступает охлаждаемая жидкость, во вторую – хладагент.
Для обеспечения однородного распределения температуры по объёму в камерах установлены мешалки. Плотность охлаждаемой жидкости 850 кг/м3, а хладагента 920 кг/м3. Объёмы камер равны и составляют 2,5 м3 каждая. Объёмный расход теплоносителя 4,12·10–3 м3/с, хладагента 5,43·10–3 м3/с.
Теплоёмкости жидкости и хладагента соответственно 3,75·103 Дж /(кг°С) и 3,14·103 Дж/(кг °С).
Поверхность теплообмена составляет 4 м2, а коэффициент теплопередачи равен K = 4360 Вт/(м2 °С). Температура охлаждаемой жидкости на входе меняется скачкообразно от 115 до 200 °С, а температура хладагента от 10 до 15 °С.
Изменение температуры по времени для каждого потока определяется решением системы:
где T10 и T20 – температуры жидкости и хладагента на входе в теплообменник.
Система уравнений (19) и (20) решается при начальных условиях:
T1 (0) = T1* ; T2 (0) = T2*.
Для определения T1* и T2*, соответствующих номинальному статическому режиму, составляются уравнения теплового баланса стационарного режима работы теплообменника:
После подстановки численных значений получаем следующую систему:
которая может быть решена численно или аналитически относительно любых двух переменных, входящих в эти уравнения.
Температуры теплоносителей на выходе теплообменника в установившемся состоянии составили: T1* = 74,46 °С, T2* = 43,93 °С.
Таким образом, получим следующую систему:
Начальные условия: Т1(0) = 74,46 °С; Т2(0) = 43,93 °С.
Графики решения системы уравнений математического описания динамики теплообменника представлены на рис. 9. На нём изображены изменения температур во времени для обоих теплоносителей на выходе теплообменника.
Температуры теплоносителей на выходе теплообменника установились через 2000 с и составили T1 = 126,6 °С; Т2 = 71,4 °С.
Численное интегрирование уравнений (21) и (22) аналогично предыдущим примерам 1 – 3. Отличие состоит в том, что интегрирование уравнений с помощью функции ode45 производится по времени.
Рис. 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));
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'); % Построение
ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ
1. Составить математическое описание теплообменного аппарата (варианты теплообменных аппаратов и исходные данные приведены в табл. 6, 7, 8).2. В зависимости от варианта смоделировать на ЭВМ статический и (или) динамический режимы теплообменника и определить температурные зависимости для всех теплоносителей.
КОНТРОЛЬНЫЕ ВОПРОСЫ
1. Чем определяется интенсивность перехода тепла от одного теплоносителя к другому?2. В чём отличие объектов идеального смешения и идеального вытеснения?
3. Что представляет собой уравнение теплового баланса (в статике и динамике) для объектов идеального смешения?
4. Что представляет собой уравнение теплового баланса (в статике и динамике) для объектов идеального вытеснения?
5. Какие функции MATLAB используются для численного интегрирования обыкновенных дифференциальных уравнений?
6. Какие дополнительные параметры содержит функция ode45 и другие аналогичные?
2,08·10–5 4,12·10– 3,44·10–5 5,21·10– 2,92·10–4 2,45·10– 5,22·10–4 2,74·10– 5,48·10–4 6,18·10– 3,12·10–5 4,05·10– 2,08·10–5 4,12·10– 3,44·10–5 5,21·10– 2,92·10–4 2,45·10– 5,22·10–4 2,74·10– 5,48·10–4 6,18·10– 3,12·10–5 4,05·10–
ПОСТРОЕНИЕ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ
АНАЛИТИЧЕСКИМ МЕТОДОМ
(МОДЕЛИРОВАНИЕ ХИМИЧЕСКИХ РЕАКЦИЙ)
Цель работы: освоить методику математического моделирования химических реакций на ЭВМ.Задания работы:
1. Ознакомиться с основными принципами моделирования химических реакций.
2. Выполнить вариант задания.
Содержание отчёта:
1. Вариант задания.
2. Математическая модель химической реакции.
3. Все расчётные формулы, используемые в работе.
4. Графики рассчитанных зависимостей для компонентов реакции.
5. Тексты программ.
КРАТКИЕ ТЕОРЕТИЧЕСКИЕ ПОЛОЖЕНИЯ
Скорость химической реакции есть изменение (уменьшение или увеличение) числа молей реагентов в результате химического взаимодействия в единицу времени на единицу объёма (для гомогенных реакций) или на единицу поверхности для гетерогенных процессов.В соответствии с этим определением скорость гомогенных химических реакции где V – объём реактора; dN/d – изменение числа молей реагентов в единицу времени.
Учитывая, что число молей N = СV, где С – концентрация реагента, в общем случае получим:
или Для реакций, идущих при постоянном объёме в замкнутой системе, второе слагаемое уравнения (25) равно нулю и уравнение приводится к виду:
Знак (+) указывает, что в реакции накапливается вещество, знак (–) – что концентрация вещества снижается.
Для гетерогенных химических реакций скорость процесса можно представить следующим соотношением:
В дальнейшем будем полагать, что в системе осуществляется гомогенный химический процесс.
Скорость химической реакции является функцией концентраций реагирующих веществ, а также температуры и иногда давления в системе:
где C – вектор концентраций реагирующих веществ; Р и Т – соответственно давление и температура.
Химические реакции можно условно разделить на простые и сложные. Химическая кинетика определяет простые реакции как такие, которые по существу проходят в одну стадию; простая реакция содержит один акт химического превращения. Скорость простой реакции зависит от концентраций только исходных веществ. В случае сложных реакций их скорость зависит от концентраций исходных, а также промежуточных (образующихся) веществ, а механизм учитывает несколько путей химического превращения (обратимые, последовательные, параллельные, разветвлённые, последовательно-параллельные реакции).
Уравнения кинетики простых и сложных химических реакций имеют некоторые различия, поэтому будем рассматривать их отдельно.
Если в простой реакции участвуют n веществ, то материальный баланс её выражается стехиометрическим уравнением:
где i – стехиометрические коэффициенты (отрицательные для исходных веществ и положительные для образующихся веществ); Ai – участвующие в реакции вещества.
Коэффициент i показывает, сколько молекул вещества Ai участвует в реакции.
Определим скорость простой химической реакции как скорость образования любого вещества, участвующего в реакции, делённую на стехиометрический коэффициент указанного компонента. При этом для исходных реагирующих компонентов стехиометрические коэффициенты принимаются со знаком минус, а для продуктов химического превращения – со знаком плюс. Так, для химической реакции (29) скорость определяется как Во многих формально-кинетических расчётах используются значения скорости реакции, выраженные через конкретные вещества. Для этих значений из формулы (30) получаем:
При использовании формулы (31) следует учитывать, что i < 0 для исходных веществ и i > 0 для продуктов реакции.
Скорость простой реакции, а также стадий сложных реакций, согласно закону действующих масс, пропорциональна наличной концентрации реагирующих веществ в некоторых степенях. Так, для реакции соответствующая зависимость будет выражаться кинетическим уравнением:
Здесь n1 – порядок реакции по веществу А, n2 – порядок по веществу В, сумму n1 + n2 называют общим, или суммарным порядком. Коэффициент пропорциональности k именуют константой скорости реакции.
Согласно уравнению Аррениуса константа скорости химической реакции зависит от температуры:
где Е – энергия активации; R – универсальная газовая постоянная; k0 – предэкспоненциальный множитель.