Подякувати Студентському архіву довільною сумою
Admin
26.02.2023 12:38
Дякуємо, що користуєтесь нашим архівом!
Інформація про навчальний заклад
ВУЗ:
Національний університет Львівська політехніка
Інститут:
Не вказано
Факультет:
Не вказано
Кафедра:
Не вказано
Інформація про роботу
Рік:
2024
Тип роботи:
Курсова робота
Предмет:
Математичне моделювання на ЕОМ
Група:
АВ
Частина тексту файла (без зображень, графіків і формул):
МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ „ЛЬВІВСЬКА ПОЛІТЕХНІКА”
КУРСОВА РОБОТА
з курсу
Математичне моделювання на ЕОМ
Варіант №15
ЛЬВІВ – 20__
ЗАВДАННЯ ДО КУРСОВОЇ РОБОТИ
1. Побудова математичної моделі та її числове дослідження
1.1. Згідно із завданням побудувати математичну модель об'єкту. Вказати всі закони, теоретичні та емпіричні залежності, які покладені в основу побудови моделі. Вказати вхідні, керуючі, збурюючі та вихідні величини, а також параметри стану системи.
1.2. Побудувати структурну схему моделі.
1.3. Для заданих у варіанті значень конструктивних параметрів, вхідних та керуючих величин числовим методом розв'язати систему відносно її параметрів стану.
1.4. Побудувати графіки розв'язків (перехідні процеси в системі). Систему розв’язати методом Рунге-Кутта з використанням зовнішніх функцій ODE23 або ODE45 пакету MATLAB.
Умовні символьні позначення:
ΔP - втрати тиску на регулюючому органі та в трубопроводі довжиною L;
P – Гідродинамічний тиск, Па;
Q – Об’ємна витрата рідини, м3/с;
T – Температура, град;
L, r – Довжина та радіус трубопроводу, м;
d – Діаметр ємності, м;
kB – Максимальна пропускна здатність регулюючого органу, м2;
ν – Кінематична в’язкість, м2/с ;
ρ – Густина рідини, Дж/кг*К ;
l – Переміщення регулюючого органу, l=[ 0 , 1 ] ;
ξ – Коефіцієнт опору рухові в турбулентному трубопроводі;
x0 – Початкова умова для параметру стану системи x ;
- Відхилення вхідної величини, керування чи параметру стану системи від номінального значення.
Умовні графічні позначення:
2. Дослідження систем шляхом лінеаризації
2.1. Показати суть методу лінеаризації в дослідженні нелінійних систем. Лінеаризувати нелінійну систему відносно одержаного стану рівноваги.
2.2. Побудувати структурну схему лінеаризованої системи та порівняти її із структурною схемою нелінійної системи, одержаною в п. 1.2.
2.3. Для заданих відхилень вхідних та керуючих величин від стану рівноваги розв’язати лінеаризовану систему, використовуючи зовнішні функції STEP. Побудувати перехідні процеси в лінеаризованій системі.
2.4. Порівняти графіки з перехідними процесами в нелінійній та лінеаризованій системах, одержаними в п. 1.2. Для цього графіки перехідних процесів в лінеаризованій та нелінійній системах для кожного параметру стану системи накласти.
3. Класичні методи дослідження систем
3.1. Лінеаризовану систему рівнянь привести до одного рівняння відносно параметрів стану:
а) Рівня в ємності h(t) ;
б) Витрати рідини в трубопроводі Q1(t);
в) Температури рідини в ємності T(t) .
Для приведення системи рівнянь до одного рівняння вищого порядку використати формулу Крамера. Записати функції передач одержаних систем.
3.2. Одержати аналітичні вирази для перехідних та імпульсних перехідних функцій систем, отриманих в п. 3.1. Для визначення коренів характеристичних рівнянь використати зовнішню функцію ROOTS.
1. Застосовуючи означення імпульсної перехідної функції та апарат узагальнених функцій.
2. Використовуючи апарат перетворення Лапласа.
Показати фізичну суть дії одиничних імпульсних вхідних сигналів на одержані лінеаризовані системи.
3.3. Побудувати графіки одержаних аналітично перехідних та імпульсних перехідних функцій систем.
3.4. Побудувати графіки перехідних та імпульсних перехідних функцій одержаних систем з використанням зовнішніх функцій STEP та IMPULSE відповідно. Одержані результати порівняти з результатами п. 3.3. Порівняння здійснити графічно шляхом накладання відповідних графіків.
3.5. Використовуючи інтеграл згортки, записати аналітичні вирази для визначення реакції системи на вхідний сигнал системи . Побудувати графіки перехідних процесів в системі
3.6. Побудувати графіки перехідних процесів в системі з використанням зовнішньої функції LSIM. Порівняти результати з одержаними в п. 3.5.
4. Частотні методи аналізу систем
4.1. Одержати аналітичні вирази для реакції системи на вхідний сигнал з врахуванням того, що . Порівняти з аналітичними виразами, одержаними в п. 3.5.
4.2. Одержати аналітичні вирази для амплітудно-частотних та фазочастотних характеристик системи.
4.3. Побудувати амплітудно-частотні та фазочастотні характеристики системи, використовуючи зовнішню функцію BODE. Порівняти графічно одержані результати з результатами п. 4.2.
5. Дослідження моделі в середовищі SIMULINK
5.1. За одержаною в п. 1.1. структурною схемою об’єкту дослідження побудувати модель у вікні моделі середовища SIMULINK. Задати згідно із завданням, збурення на вхід моделі, і одержати перехідні процеси зміни вихідних величин.
5.2. Порівняти одержані перехідні процеси з перехідними процесами з п.1.4.
Значення конструктивних параметрів:
P1=15 кПа, P2=9 кПа, P3=0,2кПа, Т1=300 К, Т2=350 К, L1=95 м, L2=6 м,
r1=0.08 м, r2=0.06 м, Кв=0.009 м2, l1=0.3, l3=0.6, ρ=1000 кг/м3, =10-5 м2/c, =0.9
Значення стрибкоподібно зміненої вхідної величини:l1=0,2*l1
1. Побудова математичної моделі та числове її дослідження
Запишемо закон збереження маси для ємності:
,
де,кг – маса рідини в першій ємності,витрати рідин.
Після підстановки m1 в (1) одержимо:
де- площа поперечного перерізу в ємності.
Вираз для невідомих значень та одержимо з емпіричних залежностей для руху рідини через відповідні регулюючі органи:
Математична модель ідеалізованого довгого ламінарного трубопроводу з регулюючим органом:
Тл
Тл=; kл=;
Для того, щоб отримати вираз для зміни температури в ємності, запишемо для цього закон збереження тепла:
,
де,Дж - тепло акумульоване в ємності,
, Дж/с – приріст тепла в першій ємності за одиницю часу ,
,Вт – потік тепла, який поступає в ємність чи виходить з неї з рідиною.
Підставимо в рівняння (6) вираз для маси рідини в ємності:
,
або ,
підставивши в ліву частину рівняння вираз (2), одержимо:
.
Об’єднавши рівняння одержимо систему рівнянь:
яка є математичною моделлю об’єкта .
1.2. Структурна схема нелінійної моделі об’єкту.
PP
h
P1
l1 Q1
Q1
h
P2
l2 T
T1
Q1 T2
h
T
1.3. Обчислення початкових умов моделі.
Початкові умови ми можемо знайти числовими методами за допомогою середовища MATLAB. Для цього створимо файл dani.m і скористаємося функцією FSOLVE.
%dani
p1=15000;
g=9.8;
p2=9000;
p3=200;
T1=300;
T2=350;
L1=95;
L2=6;
r1=0.08;
r2=0.06;
Kv=0.009;
l1=0.3;
l3=0.6;
d=0.5;
ro=1000;
ny=1e-5;
dz=0.9;
S=pi*d^2/4;
kl=(pi*r1^4)/(8*L1*ny);
kt=sqrt((pi^2*d^5)/(8*L2*dz));
tl=(L1*kl)/(pi*r1^2);
%xero
function y=xyz(x)
h=x(1);
Q1=x(2);
T=x(3);
dani;
y=[Q1+kt*sqrt((p2-ro*g*h)/ro)-Kv*l3*sqrt((ro*g*h-p3)/ro);
kl^2*((p1-ro*g*h)/ro)-Q1-(kl^2*Q1^2)/(Kv^2*l1^2);
Q1*(T1-T)+kt*sqrt((p2-ro*g*h)/ro)*(T2-T)];
%osn
x0=[0.1 ;0.1; 325];
y=fsolve('xero',x0)
Запустивши на виконання файл osn.m одержимо матрицю значень параметрів стану, серед яких виберемо останні, що будуть початковими їх значеннями:
0.00915416779408 0.00001624563996 3.44920902397511
1.4.Знаходження реакції нелінійної моделі на стрибкоподібне збурення
Після отримання початкових умов, побудуємо перехідні характеристика нашої математичної моделі. В цьому нам допоможе функція ODE23
%nelmod
function y=nelmod(t,x)
h=x(1);
Q1=x(2);
T=x(3);
dani;
l1=l1+0.2*l1;
y=[(Q1+kt*sqrt(p2/ro-g*h)-Kv*l3*sqrt(g*h-p3/ro))/S;
(kl^2*(p1/ro-g*h)-Q1-(kl*Q1^2)/(Kv^2*l1^2))/tl;
(Q1*(T1-T)+kt*sqrt(p2/ro-g*h)*(T2-T))/(S*h)];
x0=[ 0.915416779408; 0.001624563996;344.920902397507];
tol=odeset('reltol',1e-12);
[t,y]=ode23('nelmod',[0,300],x0,tol)
plot(t,y(:,1)); ylabel('h,m'); xlabel('t,c'); grid;
pause
plot(t,y(:,2)); ylabel('Q1, m^3/c'); xlabel('t,c'); grid;
pause
plot(t,y(:,3)); ylabel('T, K'); xlabel('t,c'); grid;
Результати виконання програми:
Графік перехідного процесу рівня рідини в ємності, після стрибкоподібної зміни відкриття регулюючого органу l2 від 0.4 до 0.8
Графік перехідного процесу витрати рідини в трубопроводі, після стрибкоподібної зміни відкриття регулюючого органу l2 від 0.4 до 0.8
Графік перехідного процесу температури рідини, після стрибкоподібної зміни відкриття регулюючого органу l2 від 0.4 до 0.8
2.1. Лінеаризація математичної моделі.
Початковий стан рівноваги розглядають як номінальний режим роботи об’єкту. Лінеаризують математичну модель відносно прийнятого номінального режиму.
Загальний вигляд лінеаризованої системи :
(9)
Де aij та bij - часткові похідні, які визначаються за формулами:
a11=-kt/2/S*sqrt(ro/(p2-ro*g*h0))*g-Kv*l3/2/S*sqrt(ro/(ro*g*h0-p3))*g;
a12=1/S;
a13=0;
a21=-kl^2*g/tl;
a22=-1/tl-2*kl^2*Q10/(Kv^2*l1^2*tl);
a23=0;
a31=-Q10*(T1-T0)/S/h0^2-kt*(T2-T0)/S/h0^2/sqrt(ro)*(ro*g*h0/2/sqrt(p2-ro*g*h0)+sqrt(p2-ro*g*h0))
a32=(T1-T0)/(S*h0);
a33=-Q10/S/h0-kt/S/h0*sqrt((p2-ro*g*h0)/ro);
b2=2*kl^2*Q10^2/(Kv^2*tl*l1^3);
Щоб використати для дослідження лінійної моделі матричні методи, на яких базуються функції STEP та IMPULSE, потрібно записати її в матричній формі :
або застосувавши умовні позначення матриць,
де - власна матриця системи,
- вектор параметрів стану,
- матриця вхідних величин.
Згідно із стандартною формою представлення лінійних систем, у випадку, коли рівень рідини в ємності, витрата рідини в трубопроводі та температура є вихідними величинами об’єкту, повною моделлю лінеаризованого об’єкту буде система рівнянь:
де ; .
2.2. Структурна схема лінеаризованої моделі об’єкту.
2.3. Числове дослідження реакції лінеаризованої моделі на стрибкоподібне збурення.
Реакцію лінійної моделі на одиничний стрибкоподібний вхідний сигнал визначимо за допомогою функції STEP. Для знаходження реакції на стрибкоподібну зміну переміщення
регулюючого органу від номінального значення , перехідну функцію системи, одержану з допомогою функції STEP, потрібно помножити на .
Послідовність команд для знаходження значень матриці стану А та вектора вхідних величин В, які є вхідними аргументами для функції STEP, записані в файлі matr.m та команди для побудови перехідних процесів , та , записані у файлі grafrez.m.
%matr
dani;
h0=0.915416779408;
Q10=0.001624563996;
T0=344.920902397507;
x0=[h0 Q10 T0];
tol=odeset('reltol',1e-12);
[t,y]=ode23('nelmod',[0,300],x0,tol);
a11=-kt/2/S*sqrt(ro/(p2-ro*g*h0))*g-Kv*l3/2/S*sqrt(ro/(ro*g*h0-p3))*g;
a12=1/S;
a21=-kl^2*g/tl;
a22=-1/tl-2*kl^2*Q10/(Kv^2*l1^2*tl);
a31=-Q10*(T1-T0)/S/h0^2-kt*(T2-T0)/S/h0^2/sqrt(ro)*(ro*g*h0/2/sqrt(p2-ro*g*h0)+
sqrt(p2- ro*g*h0))
a32=(T1-T0)/(S*h0);
a33=-Q10/S/h0-kt/S/h0*sqrt((p2-ro*g*h0)/ro);
b2=2*kl^2*Q10^2/(Kv^2*tl*l1^3);
A=[a11 a12 0;a21 a22 0;a31 a32 a33]
B=[0; b2;0]
C=[1 0 0;0 1 0;0 0 1]
D=[0;0;0]
Отримані результати:
A = -12.445836291144 5.092958178941 0
-0.000035117901 -0.014097135493 0
-68.802001634221 -249.919252538423 -0.088975736562
B = 0
0.0000086488293975308
0
C = 1 0 0
0 1 0
0 0 1
D = 0
0
0
%ravlin
matr;
t=[0:300];
[y,x]=step(A,B,C,D,1,t);
l10=0.3; l1=0.2*l10; lx=l1-l10;
x=lx*x
plot(t,x(:,1));grid;xlabel('t,sec');ylabel('h,m');
title('Grafik perehidnogo procesy linearuzovanoi sustemu rivnja h');
pause;
plot(t,x(:,2));grid;xlabel('t,sec');ylabel('Q1,m.kub/sec');
title('Grafik perehidnogo procesy linearuzovanoi sustemu vutratu Q1');
pause
plot(t,x(:,3));grid;xlabel('t,sec');ylabel('T,K');
title('Grafik perehidnogo procesy linearuzovanoi sustemu temperatyru T');
Графік перехідного процесу лінеаризованої системи рівня рідини в ємності.
Графік перехідного процесу лінеаризованої системи витрати рідини в трубопроводі
Графік перехідного процесу лінеаризованої системи температури рідини.
2.4. Порівняння графіків перехідних процесів, отриманих для нелінійної та лінеаризованої моделей.
%grafrez1
x0=[0.915416779408 0.001624563996 344.920902397507];
tol=odeset('reltol',1e-12);
[t1,y1]=ode45('nelmod',[0 600],x0,tol);
a=[ -12.445836291144 5.092958178941 0
-0.000035117901 -0.014097135493 0
-68.802001634221 -249.919252538423 -0.088975736562];
b= [ 0;
0.86488293975308e-5;
0];
c=[ 1 0 0;
0 1 0;
0 0 1];
d=[ 0; 0; 0];
t2=[0:600];
l10=0.3; lx=0.2*l10-l1;
[y2,x2]=step(a,b,c,d,1,t2)
x2=lx*x2
plot(t1,y1(:,1),'--',t2,x2(:,1)+0.915416779408); grid;
ylabel('h,m'); xlabel('time,sec'); pause
plot(t1,y1(:,2),'--',t2,x2(:,2)+0.001624563996); grid;
ylabel('Q,m.kub/sec'); xlabel('time,sec');
pause
plot(t1,y1(:,3),'--',t2,x2(:,3)+344.920902397507); grid;
ylabel('T'); xlabel('time,sec');
В результаті виконання програми отримаємо графіки перехідних процесів для нелінійної та лінеаризованої моделей.
Графік порівняння перехідних процесів рівня рідини в ємності для нелінійної і лінеаризованої моделей.
Графік порівняння перехідних процесів витрати рідини в трубопроводі для нелінійної і лінеаризованої моделей.
Графік порівняння перехідних процесів температури рідини для нелінійної і лінеаризованої моделей.
3. Класичні методи дослідження систем
3.1. Зведення системи лінійних рівнянь до одного диференційного рівняння вищого порядку.
За допомогою формули Крамера систему з матричної форми можна записати в диференціальній:
– формула Крамера.
3.1.а) Відносно рівня в ємності h(t):
З формули Крамера:
Визначник власної матриці системи в операторній формі:
Визначник
Підставивши вирази визначників у формулу Крамера, одержимо:
(s2-(a11+a22)s+(a11a22 - a12a21))= -a12b2;
або ввівши позначення: A2=1;A1=-(a11+a22);A0=(a11a22 – a12a21); B0= -a12b2;
Математична модель в диференціальній формі матиме вигляд:
Функція передачі для записаної системи буде мати вигляд:
3.1.б) Відносно витрати рідини в трубопроводі Q1(t):
Формула Крамера:
;
Визначник власної матриці системи в операторній формі:
Визначник
Підставивши вирази визначників у формулу Крамера, одержимо:
або ввівши позначення:
Математична модель в диференціальній формі матиме вигляд:
Функція передачі для записаної системи буде мати вигляд:
3.1.в) Відносно температури рідини Т(t):
Формула Крамера для даного випадку матиме вигляд:
,
де ;
Визначник власної матриці системи в операторній формі:
Визначник
Ввівши позначення:
Математична модель в диференціальній формі матиме вигляд:
(14)
Функція передачі для записаної системи буде мати вигляд:
;
3.2. Знаходження аналітичного виразу імпульсної перехідної та перехідної функцій системи.
3.2.а) Відносно рівня рідини в трубопроводі h(t):
Застосовуючи означення імпульсної перехідної функції,знайдемо імпульсну перехідну функцію рівня.
Після дії імпульсного сигналу (при t>0) рівняння прийме вигляд:
позначення h2(t) – це одинична імпульсна перехідна функція.
Початкову умову змінену за рахунок дії імпульсного вхідного сигналу, знайдемо, використовуючи залежність:
Для знаходження розв’язку рівняння знайдемо корінь його характеристичного рівняння:
Отже розв’язок системи матиме вигляд:
З початкових умов знаходимо коефіцієнти
Також імпульсну перехідну функцію можна знайти використовуючи перетворення Лапласа:
Отже вираз для імпульсної перехідної функції буде таким:
3.2.а) Відносно витрати рідини в ємності h(t):
Застосовуючи означення імпульсної перехідної функції знайдемо імпульсну перехідну функцію витрати рідини в ємності.
Після дії імпульсного сигналу (при t>0) рівняння прийме вигляд:
позначення h1(t) – це одинична імпульсна перехідна функція.
d1=
Для знаходження розв’язку рівняння (15) знайдемо корінь його характеристичного рівняння:
Отже розв’язок системи (15)матиме вигляд:
З початкових умов знаходимо коефіцієнти
Отже,у загальному вигляді перехідна функція буде такою:
Також імпульсну перехідну функцію можна знайти використовуючи перетворення Лапласа:
;
Отже вираз для імпульсної перехідної функції буде таким:
3.2.в) Відносно температури рідини Т(t):
Застосовуючи означення імпульсної перехідної функції знайдемо імпульсну перехідну функцію температури.
Після дії імпульсного сигналу (при t>0) рівняння прийме вигляд:
позначення h3(t)- це одинична імпульсна перехідна функція
Початкові умови змінені за рахунок дії імпульсного вхідного сигналу, знайдемо, використовуючи залежність:
Для знаходження розв’язку рівняння знайдемо корінь його характеристичного рівняння:
Отже розв’язок системи матиме вигляд:
З початкових умов знаходимо коефіцієнти
C2= ;
C3=
Також імпульсну перехідну функцію можна знайти використовуючи перетворення Лапласа:
Імпульсні перехідні та перехідні функції лінійних стаціонарних систем
Аналітичний вираз перехідної функції рівня рідини в ємності, з відповідними коефіцієнтами:
);
Аналітичний вираз перехідної функції витрати рідини в трубопроводі.
Аналітичний вираз перехідної функції температури, з відповідними коефіцієнтами.
=
3.3. Побудова графіків перехідних та імпульсних перехідних функцій, отриманих аналітично.
Побудову графіків для перехідних та імпульсних перехідних функцій, аналітичні вирази для яких отримані попередньо, виконаємо за допомогою файлу impper.m. Текст програми наведено.
%analituchnum
matr;
t=[0:1e-2:300];
%dlja rivnja
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B0=a12*b2,
k=roots([A2 A1 A0]),
C2=B0/A2/(k(2)-k(1));
C1=-C2,
h1=C1*exp(k(1)*t)+C2*exp(k(2)*t);
r1=C1/k(1)*(exp(k(1)*t)-1)+C2/k(2)*(exp(k(2)*t)-1);
%dlja vutratu
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B1=b2,B0=a11*b2,
k=roots([A2 A1 A0]),
C2=(B0*A2-A1*B1-B1*A2*k(1))/A2^2/(k(2)-k(1)),
C1=B1/A2-C2,
h2=C1*exp(k(1)*t)+C2*exp(k(2)*t);
r2=C1*(exp(k(1)*t)-1)/k(1)+C2*(exp(k(2)*t)-1)/k(2);
%dlja temperatyru
A3=1,A2=-(a11+a22+a33),A1=(a11*a22-a12*a21+a11*a33+a22*a33),A0=-a11*a22*a33+a12*a21*a33,
B1=a32*b2,B0=b2*(a12*a31-a32*a11),
k=roots([A3 A2 A1 A0]),
C3=(1/A3*(B0-A2/A3*B1)-B1*(k(2)+k(1))/A3)/(k(3)^2-k(3)*(k(2)+k(1))+k(1)*k(2)),
C2=B1/A3/(k(2)-k(1))-C3*(k(3)-k(1))/(k(2)-k(1)),
C1=-C2-C3,
h3=C1*exp(k(1)*t)+C2*exp(k(2)*t)+C3*exp(k(3)*t);
r3=C1*(exp(k(1)*t)-1)/k(1)+C2*(exp(k(2)*t)-1)/k(2)++C3*(exp(k(3)*t)-1)/k(3);
РЕЗУЛЬТАТ:
A2=1 A1 =12.45993342663699 A0 =0.17562949452371
B0 =4.404812641841658e-005
k =
-12.44582190424217
-0.01411152239482
C1 =-3.543207255112314e-006 C2 =3.543207255112314e-006
A2 =1 A1 =12.45993342663699 A0 =0.17562949452371
B1 =8.648829397530752e-006 B0 =-1.076419147916981e-004
k = -12.44582190424217
-0.01411152239482
C2 =-8.668474379743126e-006 C1 =1.731730377727388e-005
A3 =1 A2 =2.54890916319886 A1 =1.28426124867049 A0 =0.01562676363723
B1 =-0.00216150897836 B0 =-0.02993238615237
k =
-12.44582190424218
-0.08897573656186
-0.01411152239482
C3 = -0.03212871807901 C2 =0.03214844659626 C1 =-1.972851725333402e-005
3.4. Порівняння графіків перехідних, імпульсних перехідних функцій, отриманих аналітично та числовими методами.
Для перевірки правильності аналітичних розв’язків диференційних рівнянь (14) та (16) побудуємо в одній системі координат графік залежності (21),(33) та графік імпульсної перехідної функції, отриманої за допомогою функції IMPULSE, для рівня рідини в ємності h(t) та температури Т(t). Файл для порівняння графіків має ім’я por_imp.m. Текст файлу por_imp.m:
dani; matr;
t=[0:10:600];
%для рівня
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B0=a12*b2,
k=roots([A2 A1 A0]),
C2=B0/A2/(k(2)-k(1));
C1=-C2,
h1=C1*exp(k(1)*t)+C2*exp(k(2)*t);
r1=C1/k(1)*(exp(k(1)*t)-1)+C2/k(2)*(exp(k(2)*t)-1);
%для витрати
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B1=b2,B0=-a11*b2,
k=roots([A2 A1 A0]),
C2=(B0*A2-A1*B1-B1*A2*k(1))/A2^2/(k(2)-k(1)),
C1=B1/A2-C2,
h2=C1*exp(k(1)*t)+C2*exp(k(2)*t);
r2=C1*(exp(k(1)*t)-1)/k(1)+C2*(exp(k(2)*t)-1)/k(2);
%для температури
A3=1,A2=-(a11+a22+a33),A1=(a11*a22-a12*a21+a11*a33+a22*a33),A0=-a11*a22*a33+a12*a21*a33,
B1=a32*b2,B0=b2*(a12*a31-a32*a11),
k=roots([A3 A2 A1 A0]),
C3=(1/A3*(B0-A2/A3*B1)-B1*(k(2)+k(1))/A3)/(k(3)^2-k(3)*(k(2)+k(1))+k(1)*k(2)),
C2=B1/A3/(k(2)-k(1))-C3*(k(3)-k(1))/(k(2)-k(1)),
C1=-C2-C3,
h3=C1*exp(k(1)*t)+C2*exp(k(2)*t)+C3*exp(k(3)*t);
r3=C1*(exp(k(1)*t)-1)/k(1)+C2*(exp(k(2)*t)-1)/k(2)++C3*(exp(k(3)*t)-1)/k(3);
y=step(A,B,C,D,1,t);
plot(t,r1,'k',t,y(:,1),'o');grid;
xlabel('t,sec');ylabel('h,m');
title('Grafik porivnjannja perehidnux fynkcij rivnja h');pause;
plot(t,r2,'k',t,y(:,2),'o');grid;
xlabel('t,sec');ylabel('Q1,m.kub/sec');
title('Grafik porivnjannja perehidnux fynkcij vutratu Q1');pause;
plot(t,r3,'k',t,y(:,3),'o');grid;
xlabel('t,sec');ylabel('T,K');
title('Grafik porivnjannja perehidnux fynkcij temperatyru T');
Графік порівняння перехідних функцій рівня рідини в ємності знайдених аналітично та числовим методом.
Графік порівняння перехідних функцій витрати рідини
знайдених аналітично та числовим методом.
Графік порівняння перехідних функцій температури рідини знайдених аналітично та числовим методом.
Порівняння імпульсних перехідних функцій:
dani; matr;
t=[0:10:600];
%для рівня
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B0=a12*b2,
k=roots([A2 A1 A0]),
C2=B0/A2/(k(2)-k(1));
C1=-C2,
h1=C1*exp(k(1)*t)+C2*exp(k(2)*t);
r1=C1/k(1)*(exp(k(1)*t)-1)+C2/k(2)*(exp(k(2)*t)-1);
%для витрати
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B1=b2,B0=-a11*b2,
k=roots([A2 A1 A0]),
C2=(B0*A2-A1*B1-B1*A2*k(1))/A2^2/(k(2)-k(1)),
C1=B1/A2-C2,
h2=C1*exp(k(1)*t)+C2*exp(k(2)*t);
r2=C1*(exp(k(1)*t)-1)/k(1)+C2*(exp(k(2)*t)-1)/k(2);
%для температури
A3=1,A2=-(a11+a22+a33),A1=(a11*a22-a12*a21+a11*a33+a22*a33),A0=-a11*a22*a33+a12*a21*a33,
B1=a32*b2,B0=b2*(a12*a31-a32*a11),
k=roots([A3 A2 A1 A0]),
C3=(1/A3*(B0-A2/A3*B1)-B1*(k(2)+k(1))/A3)/(k(3)^2-k(3)*(k(2)+k(1))+k(1)*k(2)),
C2=B1/A3/(k(2)-k(1))-C3*(k(3)-k(1))/(k(2)-k(1)),
C1=-C2-C3,
h3=C1*exp(k(1)*t)+C2*exp(k(2)*t)+C3*exp(k(3)*t);
r3=C1*(exp(k(1)*t)-1)/k(1)+C2*(exp(k(2)*t)-1)/k(2)++C3*(exp(k(3)*t)-1)/k(3);
y=impulse(A,B,C,D,1,t);
plot(t,h1,'k',t,y(:,1),'ko');grid;
xlabel('t,sec');ylabel('h,m');title('Grafik porivnjannja impylsnux perehidnux fynkcij rivnja h');pause;
plot(t,h2,'k',t,y(:,2),'ko');grid;
xlabel('t,sec');ylabel('Q1,m.kub/sec');title('Grafik porivnjannja impylsnux perehidnux fynkcij vutratu Q1');
pause;
plot(t,h3,'k',t,y(:,3),'ko');grid;
xlabel('t,sec');ylabel('T,K');title('Grafik porivnjannja impylsnux perehidnux fynkcij temperatyru T');
В результаті виконання файлу одержано графіки:
Графік порівняння імпульсних перехідних функцій рівня рідини в ємності знайдених аналітично та числовим методом.
Графік порівняння імпульсних перехідних функцій витрати
рідини знайдених аналітично та числовим методом.
Графік порівняння імпульсних перехідних функцій температури рідини знайдених аналітично та числовим методом.
3.5. Знаходження аналітичного виразу реакції системи на вхідний сигнал .
Щоб визначити реакцію системи на вхідний сигнал , використаємо інтеграл згортки:
.
Отже за аналітичними виразами імпульсних перехідних функцій для рівня рідини в ємності та температури будемо шукати реакцію системи на заданий вхідний сигнал U(t).
для рівня рідини в ємності, з відповідними коефіцієнтами.
-
для витрати рідини, з відповідними коефіцієнтами.
для температури рідини, з відповідними коефіцієнтами.
-
3.6. Порівняння реакції системи на вхідний сигнал , отриманих за аналітичним розв’язком та числовим методом.
Реакцію системи на вхідний сигнал числовим методом можна знайти за допомогою зовнішньої функції LSIM. Для порівняння графіків перехідних процесів отриманих за аналітичним розв’язком та числовим методом скористаємося файлом sinyyys.m:
matr;
w=0.15;
t=[0:1:600];
[y]=lsim(A,B,C,D,sin(w*t),t,[0 0 0]);
%рівень
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B0=a12*b2,
k=roots([A2 A1 A0]),
C2=B0/A2/(k(2)-k(1));
C1=-C2,
N11=w*C1/(k(1)^2+w^2);
N21=w*C2/(k(2)^2+w^2);
M11=C1/sqrt(k(1)^2+w^2);
M21=C2/sqrt(k(2)^2+w^2);
fi1=atan(w/k(1));
fi2=atan(w/k(2));
y1=M11*sin(w*t+fi1)+N11*exp(k(1)*t)+M21*sin(w*t+fi2)+N21*exp(k(2)*t);
plot(t,y2,'o',t,y(:,2),'k'); grid;
xlabel('y1,m.kub/sec');ylabel('t,sec');
title('Grafik porivnjannja reakcii rivnja na sin(w`t)');pause;
%Витрата
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B1=b2,B0=a11*b2,
k=roots([A2 A1 A0]),
C2=-(B0*A2-A1*B1-B1*A2*k(1))/A2^2/(k(2)-k(1)),
C1=B1/A2-C2,
N11=w*C1/(k(1)^2+w^2);
N21=w*C2/(k(2)^2+w^2);
M11=C1/sqrt(k(1)^2+w^2);
M21=C2/sqrt(k(2)^2+w^2);
fi1=atan(w/k(1));
fi2=atan(w/k(2));
y2=M11*sin(w*t+fi1)+N11*exp(k(1)*t)+M21*sin(w*t+fi2)+N21*exp(k(2)*t);
plot(t,y1,'o',t,y(:,1),'k'); grid;
xlabel('y2,m.kub/sec');ylabel('t,sec');title('Grafik porivnjannja reakcii vutratu na sin(w`t)');pause;
%по температурі
A3=1,A2=-(a11+a22+a33),A1=(a11*a22-a12*a21+a11*a33+a22*a33),A0=-a11*a22*a33+a12*a21*a33,
B1=a32*b2,B0=b2*(a12*a31-a32*a11),
k=roots([A3 A2 A1 A0]),
C3=(1/A3*(B0-A2/A3*B1)-B1*(k(2)+k(1))/A3)/(k(3)^2-k(3)*(k(2)+k(1))+k(1)*k(2)),
C2=B1/A3/(k(2)-k(1))-C3*(k(3)-k(1))/(k(2)-k(1)),
C1=-C2-C3,
N11=w*C1/(k(1)^2+w^2);
N21=w*C2/(k(2)^2+w^2);
N31=w*C3/(k(3)^2+w^2);
M11=C1/sqrt(k(1)^2+w^2);
M21=C2/sqrt(k(2)^2+w^2);
M31=C3/sqrt(k(3)^2+w^2);
fi1=atan(w/k(1));
fi2=atan(w/k(2));
fi3=atan(w/k(3));
y3=M11*sin(w*t+fi1)+N11*exp(k(1)*t)+M21*sin(w*t+fi2)+N21*exp(k(2)*t)+M31*sin(w*t+fi3)+N31*exp(k(3)*t);
plot(t,y3,'o',t,y(:,3),'k'); grid;
xlabel('y3,K');ylabel('t,sec');title('Grafik porivnjannja reakcii temperatyru na sin(w`t)');pause;
Порівняння графіків реакції системи на вхідний сигнал для рівня рідини
Порівняння графіків реакції системи на вхідний сигнал для витрати рідини.
Порівняння графіків реакції системи на вхідний сигнал для температури рідини
4. Частотні методи аналізу систем
4.1. Визначення аналітичного виразу для реакції системи на вхідний сигнал.
Реакцію системи на заданий вхідний сигнал шукаємо у вигляді:
, (35)
Для витрати рідини в ємності Q1(t):
Позначимо: (36)
Для побудови та порівняння графіків АЧХ і ФЧХ створимо файл A4X.m в якому порівняємо графіки , знайдені аналітично і за допомогою команди BODE
matr;
A2=1,A1=-(a11+a22),A0=(a11*a22-a12*a21),B1=b2,B0=-a11*b2,
w=[0:0.1:5];
RE=(B1.*w.^2.*A1+B0.*(A0-w.^2))./((A0-w.^2).^2+(A1.*w).^2);
IM=(B1.*w.*(A0-w.^2)-A1.*B0.*w)./((A0-w.^2).^2+(A1.*w).^2);
A_W=sqrt(RE.^2+IM.^2);
[Aom,fiom]=bode([B1 B0],[A2 A1 A0],w);
plot(w,Aom,'o',w,A_W);grid;
xlabel('w.p/c'); ylabel('A(w)');pause
Fi_W=atan(IM./RE).*(180./pi);
plot(w,fiom,'o',w,Fi_W);grid;
xlabel('w.p/c'); ylabel('Fi(w)');
Графік порівняння амплітудно-частотних характеристик одержаних аналітично та за числовими методами
Графік порівняння фазо -частотних характеристик одержаних аналітично та за числовими методами
5. Дослідження моделі в середовищі SIMULINK
Встановимо параметри блоків моделі.
1. Блок стрибкоподібних вхідних сигналів Step.
В полі параметру Step time (момент прикладення стрибка сигналу) введемо
значення 0.
В полях Initial value (початкові значення) вкажемо значення вхідних величин в стані рівноваги: 0,3
В полях Final value (кінцеві значення) вкажемо значення вхідних величин після прикладення стрибка: 0,36.
2. Блоки об'єднання скалярних вхідних сигналів Mux.
В полі Number of inputs (число входів) для блоків слід ввести значення кількості входів .
3. Блоки обчислення виразу від вхідного сигналу Matlab fcn.
В полях Matlab function введемо вирази:
(u(1)+kt*sqrt(p2/ro-g*u(2))-Kv*l3*sqrt(g*u(2)-p3/ro))/S - для блоку Matlab fcn;
(kl^2*(p1/ro-g*u(1))-u(3)-(kl*u(3)^2)/(Kv^2*u(2)^2))/tl- для блоку Matlab fcn1;
(u(1)*(T1-u(3))+kt*sqrt(p2/ro-g*u(2))*(T2-u(3)))/(S*u(2))- для блоку Matlab fcn2.
4. Блоки інтегрування Integrator.
В полі Initial condition (початкове значення) введемо:
h0 – у блоці Integrator;
Q10 - у блоці Integrator1;
T0 – у блоці Integrator2.
%dani
p1=15000;
g=9.8;
p2=9000;
p3=200;
T1=300;
T2=350;
L1=95;
L2=6;
r1=0.08;
r2=0.06;
Kv=0.009;
l1=0.3;
l3=0.6;
d=0.5;
ro=1000;
ny=1e-5;
dz=0.9;
S=pi*d^2/4;
kl=(pi*r1^4)/(8*L1*ny);
kt=sqrt((pi^2*d^5)/(8*L2*dz));
tl=(L1*kl)/(pi*r1^2);
h0=0.915416779408;
Q10=0.001624563996;
T0=344.920902397507;
x0=[h0 Q10 T0];
Встановивши параметри моделювання, проміжок часу і похибку, запускаєм файл dani.m, і тільки після цього приступаємо досліджувати систему.
5.2. Порівняння графіків перехідних процесів в середовищі SIMULINK та графіків побудованих в пункті 1.4.
%Сімулінк порівняння
dani1;
tol=odeset('reltol',1e-12);
T=[0:5:300];
[t,y]=ode45('nelmod',T,x0,tol);
plot(t,y(:,1),'ko',tout,xout(:,1),'k');grid;
pause;
plot(t,y(:,2),'ko',tout,xout(:,2),'k');grid;
pause;
plot(t,y(:,3),'ko',tout,xout(:,3),'k');grid;
Результати виконання програми
Графік порівняння перехідних функцій рівня рідини в ємності знайдених аналітично та в середовищі Simulink.
Графік порівняння імпульсних перехідних функцій витрати рідини
знайдених аналітично та в середовищі Simulink.
Графік порівняння перехідних функцій температури рідини знайдених аналітично та в середовищі Simulink.
6. Висновки
На цій курсовій роботі ми пройшли весь основний матеріал викладу математичного моделювання. В цьому семестрі ми навчились досліджувати гідравлічні одно і багато ємності об’єкти, будувати перехідну, імпульсну перехідну функції та амплітудно- і фазочастотні характеристики.
Ділись своїми роботами та отримуй миттєві бонуси!
0%
Оголошення від адміністратора