Подякувати Студентському архіву довільною сумою
Admin
26.02.2023 12:38
Дякуємо, що користуєтесь нашим архівом!
Моделювання об’єктів ідеального змішування.
Інформація про навчальний заклад
ВУЗ:
Національний університет Львівська політехніка
Інститут:
Не вказано
Факультет:
Не вказано
Кафедра:
Кафедра АТХП
Інформація про роботу
Рік:
2005
Тип роботи:
Курсова робота
Предмет:
Математичне моделювання на ЕОМ
Частина тексту файла (без зображень, графіків і формул):
Міністерство освіти і науки України
Національний університет “Львівська політехніка”
Кафедра АТХП
Курсова робота
З курсу “Математичне моделювання на ЕОМ”
На тему:
“Моделювання об’єктів ідеального змішування”
Виконав:
Львів – 2005
Завдання на курсову роботу
1. Побудова математичної моделі та числове її дослідження.
1.1 Згідно із завданням побудувати математичну модель об'єкту. Вказати всі закони, теоретичні та емпіричні залежності, які покладені в основу побудови моделі. Вказати вхідні, керуючі, збурюючі та вихідні величини, а також параметри стану системи. Побудувати структурну схему моделі.
При побудові математичних моделей прийняти :
1. Всі ємності, трубопроводи та регулюючі органи теплоізольовані ;
2. Всі елементарні об'єкти є об'єктами ідеального змішування:
3. Тепло - та масообмін на границі розділу фаз рідина - повітря відсутні.
4. Фізичні властивості рідини в заданих діапазонах змін температур вважати величинами сталими . Для всіх варіантів завдань прийняти : =1000 кг/м3 ; =10-5 м2/c ; EMBED Equation.3 =.9.
5. Для математичного опису трубопроводів з ламінарною течією рідини прийняти ідеалізоване рівняння Пуазейля :
6. Для моделювання турбулентних трубопроводів використати ідеалізоване рівняння. Дарсі -Вейсбаха : EMBED Equation.2
7. У всіх варіантах завдань регулюючі органи нормально закриті з витратними характеристиками : EMBED Equation.2
8. При моделюванні довгих трубопроводів врахувати динаміку руху рідини в трубопроводах.
УМОВНІ СИМВОЛІЧНІ ПОЗНАЧЕННЯ:
ΔP - втрати тиску на регулюючому органі та в трубопроводі довжиною L;
P – Гідродинамічний тиск, Па;
Q – Об’ємна витрата рідини, м3/с;
T – Температура, град;
L, r – Довжина та радіус трубопроводу, м;
d – Діаметр ємності, м;
kB – Максимальна пропускна здатність регулюючого органу , м2 ;
ν – Кінематична в’язкість , м2/с ;
η – Динамічна в’язкість , Па*с ;
ρ – Густина рідини , Дж/кг*К ;
l – Переміщення регулюючого органу, l=[ 0 , 1 ] ;
ξ – Коефіцієнт опору рухові в турбулентному трубопроводі;
EMBED Equation.3 EMBED Equation.3 – Початкова умова для параметру стану системи x ;
EMBED Equation.3 - Відхилення вхідної величини , керування чи параметру стану стану системи від номінального значення .
УМОВНІ ГРАФІЧНІ ПОЗНАЧЕННЯ :
EMBED PBrush \* LOWER \s
EMBED PBrush \* LOWER \s
EMBED PBrush \* LOWER \s
1.2 Для заданих у варіанті значень конструктивних параметрів, вхідних та керуючих величин числовим методом розв'язати систему відносно її параметрів стану та побудувати графіки розв'язків (перехідні процеси в системі ). Систему розв’язати методом Рунге-Кутта з використанням зовнішніх функцій ODE23 та ODE45 пакету MATLAB.
2. Дослідження систем шляхом лінеаризації.
2.1 Показати суть методу лінеаризації в дослідженні нелінійних систем.
2.2 Для заданих вхідних величин та керування записати аналітичний вираз ( систему нелінійних алгебраїчних рівнянь ) для визначення параметрів стану системи в стані рівноваги. Розв'язати одержану систему нелінійних алгебраїчних рівнянь з допомогою зовнішньої функції FSOLVE пакету MATLAB.
2.3 Лінеаризувати нелінійну систему відносно одержаного стану рівноваги. Побудувати структурну схему лінеаризованої системи та порівняти її із структурною схемою нелінійної системи, одержаною в п. 1.1.
2.4 Для заданих відхилень вхідних та керуючих величин від стану рівноваги розв’язати лінеаризовану систему, використовуючи зовнішні функції ODE та STEP. Побудувати перехідні процеси в лінеаризованій системі та порівняти їх з перехідними процесами в нелінійній системі, одержаними в п. 1.2. Для порівняння графіки перехідних процесів в лінійній та нелінійній системах для кожного параметру стану системи накласти.
3. Класичні методи дослідження систем.
3.1 Лінеаризовану систему рівнянь привести до одного рівняння відносно параметрів стану :
а) Рівня в ємності h(t) ;
в) Температури рідини в ємності T(t) .
Для приведення системи рівнянь до одного рівняння вищого порядку використати формулу Крамера
Записати функції передач одержаних систем.
ЗАУВАЖЕННЯ :для варіантів завдань з двома гідравлічними ємностями систему привести до одного рівняння відносно рівня та температури в другій ємності відповідно.
3.2 Одержати аналітичні вирази для перехідних та імпульсних перехідних функцій систем, отриманих в п. 3.1 . Для визначення коренів характеристичних рівнянь використати зовнішню функцію ROOTS. Аналітичні вирази для імпульсних перехідних та перехідних функцій систем отримати двома шляхами:
1. Застосовуючи означення імпульсної перехідної функції та апарат узагальнених функцій.
2. Використовуючи апарат перетворення Лапласа.
Показати фізичну суть дії одиничних імпульсних вхідних сигналів на одержані лінеаризовані системи.
3.3 Побудувати графіки одержаних аналітично перехідних та імпульсних перехідних функцій систем.
3.4 Побудувати графіки перехідних та імпульсних перехідних функцій одержаних систем з використанням зовнішніх функцій STEP та IMPULSE відповідно. Одержані результати порівняти з результатами п. 3.3. Порівняння здійснити графічно шляхом накладання відповідних графіків.
3.5 Використовуючи інтеграл згортки записати аналітичні вирази для визначення реакції системи на вхідний сигнал системи EMBED Equation.2. Побудувати графіки перехідних процесів в системі
3.6 Побудувати графіки перехідних процесів в системі з використанням зовнішньої функції LSIM. Порівняти результати з одержаними в п. 3.5.
4. Частотні методи аналізу систем.
4.1 Одержати аналітичні вирази для реакції системи на вхідний сигнал EMBED Equation.2 з врахуванням того , що EMBED Equation.2. Порівняти з аналітичними виразами, одержаними в п. 3.5
4.2 Одержати аналітичні вирази для амплітудно-частотних та фазо-частотних характеристик системи та побудувати їх графіки.
4.3 Побудувати амплітудно-частотні та фазо-частотні характеристики системи, використовуючи зовнішню функцію BODE . Порівняти графічно одержані результати з результатами п. 4.2 .
Об’єкт моделювання:
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
Конструктивні параметри: EMBED Equation.3
Вхідні величини: EMBED Equation.3
Керуючі величини: EMBED Equation.3
Відхилення вхідних та керуючих величин: EMBED Equation.3 .
Параметри стану об’єкту: EMBED Equation.3
Вихідні величини: EMBED Equation.3
Побудова математичної моделі та числове її дослідження.
Згідно із завданням побудувати математичну модель об'єкту. Вказати всі закони, теоретичні та емпіричні залежності, які покладені в основу побудови моделі. Вказати вхідні, керуючі, збурюючі та вихідні величини, а також параметри стану системи. Побудувати структурну схему моделі.
В побудову математичної моделі покладено “Закон збереження маси” і “Закон збереження тепла”.
Згідно з законом збереження маси складаємо, диференційне рівняння зміни рівня рідини в ємності:
EMBED Equation.3 - зміна рівня рідини в ємності (1)
EMBED Equation.3 - витрата рідини на вході в об’єкту зліва; (2)
EMBED Equation.3 - витрата рідини на вході в об’єкту справа; (3)
де: EMBED Equation.3 - витратний коефіцієнт; (4)
EMBED Equation.3 - зміна витрати в довгому трубопроводі (5)
де EMBED Equation.3 (6)
EMBED Equation.3 ;- відповідні коефіцієнти (7)
EMBED Equation.3 - площа поперечного перерізу в ємності (8)
Згідно закону збереження тепла, складаємо диференціальні рівняння зміни температури рідини в ємності:
EMBED Equation.3 -зміна температури рідини в ємності (9)
На основі диференціальних рівнянь (1), (5) і (9), складаємо математичну модель даного об’єкту:
EMBED Equation.3 (10)
2. Визначення невідомих параметрів стану рівноваги.
Для знаходження реакції нелінійної моделі на стрибкоподібну зміну витрати Q, потрібно знайти значення параметрів стану в початковому стані рівноваги . Для цього перепишемо систему (10) у вигляді :
EMBED Equation.3 (11)
Систему (11) для знаходження початкових параметрів, розв’яжемо за допомогою функції ODE23, для цього ми описуємо її в файлі and_f1:
function y=and_f1(t,x);
h=x(1);Q1=x(2);T=x(3);
ro=1000;g=9.8;
p1=15000;p2=11000;p3=0;
dz=0.9;k3=0.01;
r1=0.7;r2=0.07;
L1=105;L2=11;
d=0.5;nu=1e-5;
l1=0.2;l3=0.6;
T1=300;T2=350;
%---------------------------------
s=pi*d^2/4;
k1=sqrt(4*pi^2*r1^5/(L1*dz));
k2=pi*r2^4/(8*nu*L2);
B=1+(k1/k3/l1)^2;
st=pi*r1^2;
A=k1*L1/st;
Q2=k2*(p2-ro*g*h)/ro;
Q3=k3*l3*sqrt((ro*g*h-p3)/ro);
%---------------------------------
y=[(Q1+Q2-Q3)/s;
(k1^2*(p1-ro*g*h)/ro-B*Q1^2)/A;
1/s*((Q1*(T1-T)+Q2*(T2-T))/h); ];
Для запуску цієї програми в середовищі MATLAB створюємо файл and_gr.m:
TT=[0 5];
y0=[0.01 0.00001 3.39]*100;
tol=odeset('reltol',1e-12,'AbsTol',[1e-12, 1e-16, 1e-12]);
[t,y]=ode45('and_f1',TT,y0,tol);
plot(t,y(:,1));grid;pause;
plot(t,y(:,2));grid;pause;
TT=[0 80];
y0=[0.01103825609477 0.00004090120177 3.39636855544612]*100;
tol=odeset('reltol',1e-12,'AbsTol',[1e-12, 1e-16, 1e-12]);
[t,y]=ode45('and_f1',TT,y0,tol);
plot(t,y(:,3));grid;pause;close
Після запуску програми на виконання, я отримав слідуючі значення початкових параметрів стану:
» y
y = 1.0e+002 *
0.01103825609477 0.00004090120177 3.39636855544612
Графіки перехідних процесів h(t), Q1(t), Т(t) отримаємо підставивши в файл and_gr значення тиску Р2=1.2*P2
function y=and_f1(t,x);
h=x(1);Q1=x(2);T=x(3);
ro=1000;g=9.8;
p1=15000;p2=13200;p3=0;
dz=0.9;k3=0.01;
r1=0.7;r2=0.07;
L1=105;L2=11;
d=0.5;nu=1e-5;
l1=0.2;l3=0.6;
T1=300;T2=350;
%---------------------------------
s=pi*d^2/4;
k1=sqrt(4*pi^2*r1^5/(L1*dz));
k2=pi*r2^4/(8*nu*L2);
B=1+(k1/k3/l1)^2;
st=pi*r1^2;
A=k1*L1/st;
Q2=k2*(p2-ro*g*h)/ro;
Q3=k3*l3*sqrt((ro*g*h-p3)/ro);
%---------------------------------
y=[(Q1+Q2-Q3)/s;
(k1^2*(p1-ro*g*h)/ro-B*Q1^2)/A;
1/s*((Q1*(T1-T)+Q2*(T2-T))/h); ];
Запускаючий файл and_gr.m.
TT=[0 5];
y0=[0.01103825609477 0.00004090120177 3.39636855544612]*100;
tol=odeset('reltol',1e-12,'AbsTol',[1e-12, 1e-16, 1e-12]);
[t,y]=ode45('and_f1',TT,y0,tol);
plot(t,y(:,1));grid;pause;
plot(t,y(:,2));grid;pause;
TT=[0 80];
y0=[0.01103825609477 0.00004090120177 3.39636855544612]*100;
tol=odeset('reltol',1e-12,'AbsTol',[1e-12, 1e-16, 1e-12]);
[t,y]=ode45('and_f1',TT,y0,tol);
plot(t,y(:,3));grid;pause;close
Запустивши програми на виконання, я отримав наступні графіки зміни параметрів стану в залежності від часу:
Рисунок 1 Перехідний процес h(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20.
Рисунок 2 Перехідний процес Q(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20.
Рисунок 3 Перехідний процес Q(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20,
3. Дослідження систем шляхом лінеаризації.
Початковий стан рівноваги будемо розглядати як номінальний режим роботи об’єкту. Отже h0=hH ;Q0 =Qн; T0=TH ; і т. д. Ліреаризуємо математичну модель відносно прийнятого номінального режиму :
Перепишемо ще раз систему (10):
EMBED Equation.3 (12)
Тоді ліанеризована модель буде мати вигляд :
EMBED Equation.3 (13)
де:
EMBED Equation.3 ,
EMBED Equation.3
EMBED Equation.3
вихідні величини : y1(t)=h(t)
y2(t)=Q1(t)
y3(t)=T(t)
Отже лінеаризована модель буде мати вигляд:
EMBED Equation.3 де : EMBED Equation.3 ; EMBED Equation.3 ; EMBED Equation.3 ;
EMBED Equation.3 ; EMBED Equation.3 EMBED Equation.3 ; EMBED Equation.3; EMBED Equation.3 .
EMBED Equation.3 EMBED Equation.3
EMBED Equation.3 EMBED Equation.3
EMBED Equation.3 EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
Отже:
EMBED Equation.3 B= EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
Реакцію лінійної моделі на одиничний стрибкоподібний сигнал можна визначити за допомогою функції Step. Для знаходження реакції на стрибкоподібне відхилення вхідної величини від нмінального значення EMBED Equation.3 ,перехідну функцію системи, одержану за допомогою функції Step потрібно домножити на EMBED Equation.3 .
Послідовність команд для знаходження значень матриці стану А та вектора
вхідних величин b, які є вхідними аргументами функції STEP, записана у файлі and_matr:
ro=1000;g=9.8;
p1=15000;p2=11000;p3=0;
dz=0.9;k3=0.01;
r1=0.7;r2=0.07;
L1=105;L2=11;
d=0.5;nu=1e-5;
l1=0.2;l3=0.6;
T1=300;T2=350;
%---------------------------------
s=pi*d^2/4;
k1=sqrt(4*pi^2*r1^5/(L1*dz));
k2=pi*r2^4/(8*nu*L2);
B=1+(k1/k3/l1)^2;
st=pi*r1^2;
A=k1*L1/st;
h0=1.103825609477;
Q10=0.004090120177;
T0=339.636855544631;
a11=(-k2*g-1/2*k3*l3*g/((ro*g*h0-p3)/ro)^(1/2))/s;
a12=1/s;
a13=0;
a21=-k1^2*g/A;
a22=-2*B*Q10/A;
a23=0;
a31=-k2*g*(T2-T0)/h0/s-(Q10*(T1-T0)+k2*(p2-ro*g*h0)/ro*(T2-T0))/s/h0^2;
a32=1/s*(T1-T0)/h0;
a33=1/s*(-Q10-k2*(p2-ro*g*h0)/ro)/h0;
b1=k2/ro/s;
b2=0;
b3=1/s*k2/ro*(T2-T0)/h0;
a=[a11 a12 a13;
a21 a22 a23;
a31 a32 a33 ]
b=[b1; b2; b3]
c=[1 0 0;
0 1 0;
0 0 1 ]
d=[0; 0; 0]
Запустивши програму на виконання, я отримала наступні результати:
a =
1.0e+002 *
-0.04323670896405 0.05092958178941 0
-0.00038070781638 -0.07945048116012 0
-0.40164894677099 -1.82881105403202 -0.00091050883718
b =
0.00043654545455
0
0.00409845864052
c =
1 0 0
0 1 0
0 0 1
d =
0
0
0
Для порівняння перехідних процесів потрібно накласти в одній системі координат графіки нелінійної та лінійної моделі. Для цього було створюємо файл and_lin
x0=[ 0.01103825609477 0.00004090120177 3.39636855544631]*100;
p20=11000; px=0.2*p20;
TT=[0:0.1:5];
tol=odeset('RelTol',1e-12,'AbsTol',[1e-8 1e-8 1e-8]);
[t1,y1]=ode45('and_f1',TT,x0,tol);
%-----------------------
and_matr;
t=[0:0.1:5];
[y,x]=step(a,b,c,d,1,t);
x=px*x;
%------------------------
plot(t1,y1(:,1),'o',t,x(:,1)+ 1.103825609477);grid;
ylabel('h,m'); xlabel('time,sec'); pause;
plot(t1,y1(:,2),'o',t,x(:,2)+ 0.004090120177);grid;
ylabel('Q,m.kub/sec'); xlabel('time,sec');
pause;
plot(t1,y1(:,3),'o',t,x(:,3)+339.636855544631);grid;
ylabel('T,K'); xlabel('time,sec');pause;close;
Врезультаті виконання програми отримаємо графіки перехідних процесів в нелінійній та лінеарезованій моделях:
Рисунок 4 Порівняльний графік перехідного процесу h(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20, в не лінійній та лінеаризованій моделях
.
Рисунок 5 Порівняльний графік перехідного процесу Q1(t) з номінального початкового стану рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20, в не лінійній та лінеаризованій моделях
Рисунок 6 Порівняльний графік перехідного процесу T(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20, в не лінійній та лінеаризованій моделях
---- нелінійна модель
лінеаризована модель
4.Класичні методи дослідження системи
Запишимо лінеаризовану математичну модель:
EMBED Equation.3 (14)
у вигляді:
де EMBED Equation.2 ,
EMBED Equation.2 .
Початкові умови: EMBED Equation.3
Зведення системи лінійних рівнянь до одного диференціального рівняння.
Зведемо систему лінійних рівнянь до одного диференційного рівняння .В класичній формі система матиме вигляд:
EMBED Equation.3 (15)
Запишемо систему (14) в операторній формі. Для цього введемо діагональну матрицю диференціювання EMBED Equation.3 . Тоді, система (15) прийме вигляд
EMBED Equation.3 ,
EMBED Equation.3 , EMBED Equation.3 , EMBED Equation.3 .
За правилом Крамера:
EMBED Equation.3 , (16)
EMBED Equation.3 .
Визначник власної матриці системи в операторній формі:
EMBED Equation.3
Підставивши вирази визначників у формулу (2), одержимо:
EMBED Equation.3
EMBED Equation.3 .
Або, ввівши позначення EMBED Equation.3 , EMBED Equation.3 , отримаємо кінцевий вигляд її моделі в операторній формі:
EMBED Equation.3 ; (17)
в диференціальній формі
EMBED Equation.3
Знаходження імпульсної перехідної функції системи.
Для лінійної математичної моделі справедливі принципи суперпозиції, однорідності та комутативності, тому загальну реакцію моделі можна шукати як суму реакції на збурення
EMBED Equation.3 та EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
знайдемо реакцію r1(t) на стрибкоподібне збурення EMBED Equation.3
EMBED Equation.3 + EMBED Equation.3
EMBED Equation.3
заміна EMBED Equation.3 EMBED Equation.3 дозволить звести диференційне рівняння до однорідного:
EMBED Equation.3
Для даного варіанту EMBED Equation.3 47,2259; EMBED Equation.3 14,8423;
Характеристичне рівняння отримане за диф. рівнянням - EMBED Equation.3
EMBED Equation.3 , EMBED Equation.3 . EMBED Equation.3
Оскільки корені характеристичного рівняння дійсні різні, то розв’язок диф. рівнянням має вигляд:
EMBED Equation.3 , перейшовши від z до r1 маємо EMBED Equation.3 загальний розв’язок диф. рівняння має вигляд EMBED Equation.3
EMBED Equation.3 r(t)= EMBED Equation.3 ; EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3
EMBED Equation.3 ; EMBED Equation.3
Отримаємо: EMBED Equation.3 EMBED Equation.3 (18)
Невідомі коефіцієнти EMBED Equation.3 і EMBED Equation.3 EMBED Equation.3 знайдемо з початкових умов EMBED Equation.3 , EMBED Equation.3 Значення EMBED Equation.3 знаходимо з першого рівняння системи:
EMBED Equation.3
враховуючи ,що при t=(0+) ; EMBED Equation.3 , EMBED Equation.3 , отримаємо:
EMBED Equation.3
EMBED Equation.3 ;
EMBED Equation.3
при EMBED Equation.3
EMBED Equation.3
EMBED Equation.3 , EMBED Equation.3 EMBED Equation.3 , де EMBED Equation.3 . EMBED Equation.3
EMBED Equation.3 EMBED Equation.3
Порівняння графіків імпульсних перехідних функцій, отриманих за аналітичним розв’язком та числовим методом.
Для перевірки правильності аналітичного розв’язку диференційного рівняння (17) побудуємо в одній системі координат графік залежності (18) та графік перехідної функції, отриманої за допомогою функції.STEP.Послідовність команд, яку потрібно виконати для побудови графіків, зібрана у script-файлі and_klsmet:
Px=1; t=[0:3:200];
a=[ -0.04323670896405 0.05092958178941 0;
-0.00038070781638 -0.07945048116012 0;
-0.48333403299959 -1.82881105403202 -0.00961115857528 ];
b =[0.00043654545455;
0;
0.00409845864052 ];
A1=-(a(1,1)+a(2,2)); A2=a(1,1)*a(2,2)-a(1,2)*a(2,1);
A3=A1/A2; A4=1/A2;
D=sqrt(A3^2-4*A4);
lb1=(-A3+D)/(2*A4); lb2=(-A3-D)/(2*A4);
N1=b(1,1)*(lb1-a(2,2)); N2=b(1,1)*(lb2-a(2,2)); N3=a(2,2)*b(1,1)*(Px/A2);
C2=(b(1,1)*Px-lb1*N3)/(lb2-lb1)/N2;
C1=(N3-N2*C2)/N1;
r=N1*C1*exp(lb1*t)+N2*C2*exp(lb2*t)-N3;
%Розв"язування лін.с-ми числоим методом
c=[1 0 0;0 1 0;0 0 1]; d=[0;0;0];
tt=[1:1:200];
[y,x]=step(a,b,c,d,1,t);
plot(t,r,'o',t,x(:,1));grid;
title('Perexidna f');
ylabel('h-h0,m');xlabel('time,sec');pause;
Рис.12. Порівняння графіків перехідної функції:
- за аналітичною залежністю; о – із застосуванням функції STEP.
Знаходження імпульсної перехідної функції системи.
Для лінійної моделі імпульсна перехідна функція h(t)=dr(t)/dt, тобто
EMBED Equation.3 (19)
Порівняння графіків імпульсної перехідної функції.
Для порівняння накладемо графік перехідної функції, побудованої за аналітичною залежністю та графік імпульсної перехідної функції знайденої числовим методом в одній системі координат. Для числового знаходження реакції на EMBED Equation.3 використаємо функцію IMPULSE.
Для побудови графіків перехідних функцій додамо у файл and_klsmet такі стрічки:
h=lb1*N1*C1*exp(lb1*tt)+lb2*N2*C2*exp(lb2*tt);
[y,x]=impulse(a,b,c,d,1,tt);
plot(tt,h,'o',tt,x(:,1),'k');grid;
title('impulsna perexidna f');
ylabel('h-h0,m'); xlabel('time,sec');pause;close;
Після виконання модифікованого файлу and_klsmet одержимо графік:
Рис.13. Порівняння графіків перехідної функції:
- за аналітичною залежністю;
о – із застосуванням функції IMPULSE.
5. Аналітичні вирази для визначення реакції системи на вхідний синусоїдальний сигнал ( EMBED Equation.2 ) за допомогою інтеграла згортки
Реакцію системи (аналітичний її вираз) на вхідний сигнал EMBED Equation.2 знайдемо за допомогою iнтеграла згортки:
EMBED Equation.2 .
EMBED Equation.3
EMBED Equation.3
Порівняємо графіки, отримані за допомогою аналітичних виразів та з допомогою функції MatLab lsim. Відповідні програми знаходяться у файлі:
and_matr;
A2=1; A1=-(a11+a22); A0=a11*a22-a12*a21;
B1=b(1,1);B0=a(2,2)*b(1,1);
w=0.5;
t=[0:0.1:50];
A1=-(a(1,1)+a(2,2)); A2=a(1,1)*a(2,2)-a(1,2)*a(2,1);
A3=A1/A2; A4=1/A2;
D=sqrt(A3^2-4*A4);
lb1=(-A3+D)/(2*A4); lb2=(-A3-D)/(2*A4);
d0=1/A2*B1; d1=1/A2*(B0-A1*d0);
C2=(d1-d0*lb1)/(lb2-lb1);C1=d0-C2;
bs=-(C1*w./(lb1^2+w.^2)+C2*w./(lb2^2+w.^2));
as=(-lb1*C1)./(lb1^2+w.^2)+(-lb2*C2)./(lb2^2+w.^2);
K=C1*w*exp(lb1*t)/(lb1^2+w.^2)+C2*w*exp(lb2*t)/(lb2^2+w.^2);
fi=atan(bs/as);
y=(sqrt(as.^2+bs.^2).*sin(w.*t+fi)*68+K);
sys=ss(a,b,c,d);
x=lsim(sys,sin(w.*t),t);
plot(t,x(:,1),'or',t,y);grid;pause;close;
Запустивши програму на виконання, одержимо:
Рис.14. Порівняння графіків перехідних процесів:
- за аналітичною залежністю;о – із застосуванням функції LSIM
6. Частотні методи аналізу системи (амплітудо- та фазочастотні характеристики)
Амплітудочастотна характеристика (А(w)) - це залежність відношення амплітуди вихідного періодичного сигналу до вхідної амплітуди цього сигналу як функції частоти в усталеному режимі роботи.
Фазочастотна характеристика ((w)) - це залежність зміни фази вихідного сигналу по відношенню до фази вхідного сигналу від частоти (різниця фаз (w)=вих-вх) в усталеному режимі роботи.
Знайдемо амплітудочастотну та фазочастотну характеристику системи, зведеної відносно рівня в другій ємності:
EMBED Equation.2 .
EMBED Equation.3
EMBED Equation.3
Побудуємо амплітудно-частотні та фазо-частотні характеристики системи, використовуючи зовнішню функцію BODE і порівняємо графічно одержані результати з результатами одержаними за аналітичними залежностями. Для цього створюємо файл:
and_matr;
A2=1; A1=-(a11+a22); A0=a11*a22-a12*a21;
B1=b1;B0=a22*b1; w=[0:5:200];
H=(B0+B1*i*w)./(-A2*w.^2+A1*i*w+A0);
RE=real(H);
IM=imag(H);
A_W=sqrt(RE.^2+IM.^2);
plot(w,A_W);grid;pause;
H=(B0+B1*i*w)./(-A2*w.^2+A1*i*w+A0);
RE=real(H);
IM=imag(H);
Fi_W=atan(IM./RE)*180/pi;
plot(w,Fi_W);grid;pause;
[Aom,fiom]=bode([B1 B0],[A2 A1 A0],w)
plot(w,Aom,'ro',w,A_W);grid;
xlabel('w.p/c'); ylabel('A(w)');pause;
plot(w,fiom,'ro',w,Fi_W);grid;
xlabel('w.p/c');
ylabel('Fi(w)');pause;
Рис. 15 Порівняння графіків амплітудночастотної характеристики системи, зведеної до одного диференціального рівняння відносно h, отримані з допомогою: 1 - аналітичного виразу; 2 - функції MatLab bode.
Рис. 16 Порівняння графіків фазочастотної характеристики системи, зведеної до одного диференціального рівняння відносно h, отримані з допомогою:
1 - аналітичного виразу; 2 - функції MatLab bode
7.Дослідження моделі в середовищі SIMULINK
Побудова структурної схеми моделі.
Структурна схема моделі, яка відображає причинно-наслідкові зв’язки між входами та виходами системи:
SHAPE \* MERGEFORMAT EMBED Equation.3
EMBED Equation.3
EMBED Equation.3
h
l
Q
T
P2
P3
P1
Рис. 7 Структурне зображення математичної моделі об’єкту
Побудова структурної схеми моделі у вікні SIMULINK.
Копіюю блоки – елементи моделі з бібліотек у вікно системи:
Рис.8 Cтруктурна схема моделі у вікні SIMULINK.
Встановлення параметрів блоків моделі.
1.Блоки стрибкоподібних вхідних сигналів STEP
В полі параметру Step time введемо -0;
В полях Initial value вкажемо значення вхідних величин в стані рівноваги: 11000 і збурення в полі Final value 11000*1.2.
У всіх блоці Step в полі Sample time – 0;
2.Блоки об’єднання скалярних вхідних сигналів Mux
В полі Namber of inputs – 2.
3. Блоки обчислення виразу від вхідного сигналу Step Fcn.
В полях Matlab function введемо вирази:
u(2)+k2*(u(3)-ro*g*u(1))/ro-k3*l3*sqrt((ro*g*u(1)-p3)/ro)-для блоку Matlab function;
k1^2*(p1-ro*g*u(1))/ro-B*u(2)^2 - для блоку Matlab function1;
(u(1)*(T1-u(3))+(k2*(p2-ro*g*u(2))/ro)*(T2-u(3)))/u(2)-для блоку Matlab function2;
5. Блоки підсилення Gain:
1/S – Gain;
1/A – Gain1;
1/Sh0 – Gain2;
6. Блоки інтегрування Integrator.
В полі Initial condition введемо:
h0 – у блоці Integrator
Q10 – у блоці Integrator1
Т0 – у блоці Integrator2
Створимо окремий файл розрахунку коефіцієнтів моделі та початкових значень параметрів стану:
and_dani_sim
ro=1000;g=9.8;
p1=15000;p2=11000*1.2;p3=0;
dz=0.9;k3=0.01;
r1=0.7;r2=0.07;
L1=105;L2=11;
d=0.5;nu=1e-5;
l1=0.2;l3=0.6;
T1=300;T2=350;
%---------------------------------
s=pi*d^2/4;
k1=sqrt(4*pi^2*r1^5/(L1*dz));
k2=pi*r2^4/(8*nu*L2);
B=1+(k1/k3/l1)^2;
st=pi*r1^2;
A=k1*L1/st;
Встановлення параметрів моделювання.
Відкривши вікно параметрів моделювання в полях start time і stop time введемо значення початку і кінця моделювання: Т=[0 20].
5.Дослідження динаміки зміни параметрів стану, спричинені стрибкоподібними змінами вхідних величин від стану рівноваги. Перед початком моделювання потрібно визначити в робочій області Matlab значення змінних, не заданих у числовій формі в моделі. Для цього з командного рядка Matlab виконаємо команду:
>> sml
6.Порівняння перехідних процесів із результатами, одержаними шляхом розв’язування системи нелінійних диференційних рівнянь.
Для цього виконаємо у вікні Matlab команди:
>>x0=[0.00726081907601 0.00004446455541 3.39327546235733]*100;
>> [t,y]=ode45('stl_r15',[0 20],x0);
де and_f1 функція для обчислення значень вектора похідних в алгоритмі ode45
function y=and_f1(t,x);
h=x(1);Q1=x(2);T=x(3);
ro=1000;g=9.8;
p1=15000;p2=13200;p3=0;
dz=0.9;k3=0.01;
r1=0.7;r2=0.07;
L1=105;L2=11;
d=0.5;nu=1e-5;
l1=0.2;l3=0.6;
T1=300;T2=350;
%---------------------------------
s=pi*d^2/4;
k1=sqrt(4*pi^2*r1^5/(L1*dz));
k2=pi*r2^4/(8*nu*L2);
B=1+(k1/k3/l1)^2;
st=pi*r1^2;
A=k1*L1/st;
Q2=k2*(p2-ro*g*h)/ro;
Q3=k3*l3*sqrt((ro*g*h-p3)/ro);
%---------------------------------
y=[(Q1+Q2-Q3)/s;
(k1^2*(p1-ro*g*h)/ro-B*Q1^2)/A;
1/s*((Q1*(T1-T)+Q2*(T2-T))/h); ];
Для порівняння результатів побудуємо графіки. Для цього використовуємо файл and_gr_sim
y0=[0.01103825609477 0.00004090120177 3.39636855544612]*100;
t0=0; tf=10;
options = odeset('RelTol',1e-4,'AbsTol',1e-7);
[t1,y]=ode23('and_f1',[t0 tf],y0,options);
%************************************************
plot( t1 , y(:,1) ,tout , xout(:,1),'or' );grid;pause;
xlabel('t.sec');
ylabel('h.m');
plot(t1, y(:,2), tout, xout(:,2),'or');grid;pause;
xlabel('t.sec');
ylabel('Q1.m.kub/sec');
plot(t1, y(:,3), tout, xout(:,3),'or');grid;pause;
xlabel('t.sec');
ylabel('T.K');
Рис.9 Порівняльний графік перехідного процесу h(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20, в не лінійній моделі та результатами отриманими в середовищі SIMULINK.
Рис.10 Порівняльний графік перехідного процесу Q(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20, в не лінійній моделі та результатами отриманими в середовищі SIMULINK.
Рис.11 Порівняльний графік перехідного процесу T(t) з номінального початкового стану
рівноваги до нового стану рівноваги, при нанесенні на об’єкт збурення тиском Р2=1.2*Р20, в не лінійній моделі та результатами отриманими в середовищі SIMULINK.
Висновки
Виконавши курсову роботу з дисципліни «Математичне моделювання на ЕОМ» я закріпив знання, отримані на лекціях та практичних заняттях. Навчився на прикладі системи проточних ємностей досліджувати її властивості, зокрема знаходити реакцію системи на стрибкоподібне та імпульсне збурення, будувати перехідну, імпульсну перехідну функції та амплітудно- і фазочастотні характеристики .
В процесі виконання цієї курсової роботи я навчився також користуватися функціями програмного пакету MatLab, який значно прискорює проведення складних розрахунків.
Ділись своїми роботами та отримуй миттєві бонуси!
0%
Оголошення від адміністратора