Подякувати Студентському архіву довільною сумою
Admin
26.02.2023 12:38
Дякуємо, що користуєтесь нашим архівом!
Моделювання об’єктів ідеального змішування
Інформація про навчальний заклад
ВУЗ:
Національний університет Львівська політехніка
Інститут:
Не вказано
Факультет:
Не вказано
Кафедра:
Кафедра АТХП
Інформація про роботу
Рік:
2005
Тип роботи:
Курсова робота
Предмет:
Математичне моделювання на ЕОМ
Група:
К
Частина тексту файла (без зображень, графіків і формул):
Міністерство освіти і науки України
Національний університет “Львівська політехніка”
Кафедра АТХП
Курсова робота
З курсу “Математичне моделювання на ЕОМ”
На тему:
“Моделювання об’єктів ідеального змішування”
Завдання на курсову роботу
1. Побудова математичної моделі та числове її дослідження.
1.1 Згідно із завданням побудувати математичну модель об'єкту. Вказати всі закони, теоретичні та емпіричні залежності, які покладені в основу побудови моделі. Вказати вхідні, керуючі, збурюючі та вихідні величини, а також параметри стану системи. Побудувати структурну схему моделі.
При побудові математичних моделей прийняти :
1. Всі ємності, трубопроводи та регулюючі органи теплоізольовані ;
2. Всі елементарні об'єкти є об'єктами ідеального змішування:
3. Тепло - та масообмін на границі розділу фаз рідина - повітря відсутні.
4. Фізичні властивості рідини в заданих діапазонах змін температур вважати величинами сталими . Для всіх варіантів завдань прийняти : ( =1000 кг/м3 ; (=10-5 м2/c ; =.9.
5. Для математичного опису трубопроводів з ламінарною течією рідини прийняти ідеалізоване рівняння Пуазейля :
6. Для моделювання турбулентних трубопроводів використати ідеалізоване рівняння. Дарсі -Вейсбаха :
7. У всіх варіантах завдань регулюючі органи нормально закриті з витратними характеристиками :
8. При моделюванні довгих трубопроводів врахувати динаміку руху рідини в трубопроводах.
УМОВНІ СИМВОЛІЧНІ ПОЗНАЧЕННЯ:
ΔP - втрати тиску на регулюючому органі та в трубопроводі довжиною L;
P – Гідродинамічний тиск, Па;
Q – Об’ємна витрата рідини, м3/с;
T – Температура, град;
L, r – Довжина та радіус трубопроводу, м;
d – Діаметр ємності, м;
kB – Максимальна пропускна здатність регулюючого органу , м2 ;
ν – Кінематична в’язкість , м2/с ;
η – Динамічна в’язкість , Па*с ;
ρ – Густина рідини , Дж/кг*К ;
l – Переміщення регулюючого органу, l=[ 0 , 1 ] ;
ξ – Коефіцієнт опору рухові в турбулентному трубопроводі;
– Початкова умова для параметру стану системи x ;
- Відхилення вхідної величини , керування чи параметру стану стану системи від номінального значення .
УМОВНІ ГРАФІЧНІ ПОЗНАЧЕННЯ :
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 Використовуючи інтеграл згортки записати аналітичні вирази для визначення реакції системи на вхідний сигнал системи . Побудувати графіки перехідних процесів в системі
3.6 Побудувати графіки перехідних процесів в системі з використанням зовнішньої функції LSIM. Порівняти результати з одержаними в п. 3.5.
4. Частотні методи аналізу систем.
4.1 Одержати аналітичні вирази для реакції системи на вхідний сигнал з врахуванням того , що . Порівняти з аналітичними виразами, одержаними в п. 3.5
4.2 Одержати аналітичні вирази для амплітудно-частотних та фазо-частотних характеристик системи та побудувати їх графіки.
4.3 Побудувати амплітудно-частотні та фазо-частотні характеристики системи, використовуючи зовнішню функцію BODE . Порівняти графічно одержані результати з результатами п. 4.2 .
Об’єкт моделювання:
Конструктивні параметри:
Вхідні величини:
Керуючі величини:
Відхилення вхідних та керуючих величин: .
Параметри стану об’єкту:
Вихідні величини:
Побудова математичної моделі та числове її дослідження.
Згідно із завданням побудувати математичну модель об'єкту. Вказати всі закони, теоретичні та емпіричні залежності, які покладені в основу побудови моделі. Вказати вхідні, керуючі, збурюючі та вихідні величини, а також параметри стану системи. Побудувати структурну схему моделі.
В побудову математичної моделі покладено “Закон збереження маси” і “Закон збереження тепла”.
Згідно з законом збереження маси складаємо, диференційне рівняння зміни рівня рідини в ємності:
- зміна рівня рідини в ємності (1)
- витрата рідини на вході в об’єкту зліва; (2)
- витрата рідини на вході в об’єкту справа; (3)
де:- витратний коефіцієнт; (4)
- зміна витрати в довгому трубопроводі (5)
де (6)
;- відповідні коефіцієнти (7)
- площа поперечного перерізу в ємності (8)
Згідно закону збереження тепла, складаємо диференціальні рівняння зміни температури рідини в ємності:
-зміна температури рідини в ємності (9)
На основі диференціальних рівнянь (1), (5) і (9), складаємо математичну модель даного об’єкту:
(10)
2. Визначення невідомих параметрів стану рівноваги.
Для знаходження реакції нелінійної моделі на стрибкоподібну зміну витрати Q, потрібно знайти значення параметрів стану в початковому стані рівноваги . Для цього перепишемо систему (10) у вигляді :
(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):
(12)
Тоді ліанеризована модель буде мати вигляд :
(13)
де:
,
вихідні величини : y1(t)=h(t)
y2(t)=Q1(t)
y3(t)=T(t)
Отже лінеаризована модель буде мати вигляд:
де : ; ; ;
; ; ; .
Отже:
B=
Реакцію лінійної моделі на одиничний стрибкоподібний сигнал можна визначити за допомогою функції Step. Для знаходження реакції на стрибкоподібне відхилення вхідної величини від нмінального значення,перехідну функцію системи, одержану за допомогою функції Step потрібно домножити на.
Послідовність команд для знаходження значень матриці стану А та вектора
вхідних величин 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.Класичні методи дослідження системи
Запишимо лінеаризовану математичну модель:
(14)
у вигляді:
де ,
.
Початкові умови:
Зведення системи лінійних рівнянь до одного диференціального рівняння.
Зведемо систему лінійних рівнянь до одного диференційного рівняння .В класичній формі система матиме вигляд:
(15)
Запишемо систему (14) в операторній формі. Для цього введемо діагональну матрицю диференціювання . Тоді, система (15) прийме вигляд
,
, , .
За правилом Крамера:
, (16)
.
Визначник власної матриці системи в операторній формі:
Підставивши вирази визначників у формулу (2), одержимо:
.
Або, ввівши позначення , , отримаємо кінцевий вигляд її моделі в операторній формі:
; (17)
в диференціальній формі
Знаходження імпульсної перехідної функції системи.
Для лінійної математичної моделі справедливі принципи суперпозиції, однорідності та комутативності, тому загальну реакцію моделі можна шукати як суму реакції на збурення
та
знайдемо реакцію r1(t) на стрибкоподібне збурення
+
заміна дозволить звести диференційне рівняння до однорідного:
Для даного варіанту 47,2259;14,8423;
Характеристичне рівняння отримане за диф. рівнянням -
, .
Оскільки корені характеристичного рівняння дійсні різні, то розв’язок диф. рівнянням має вигляд:
, перейшовши від z до r1 маємо загальний розв’язок диф. рівняння має вигляд
r(t)=;
;
Отримаємо: (18)
Невідомі коефіцієнти і знайдемо з початкових умов , Значення знаходимо з першого рівняння системи:
враховуючи ,що при t=(0+) ; ,, отримаємо:
;
при
, , де .
Порівняння графіків імпульсних перехідних функцій, отриманих за аналітичним розв’язком та числовим методом.
Для перевірки правильності аналітичного розв’язку диференційного рівняння (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, тобто
(19)
Порівняння графіків імпульсної перехідної функції.
Для порівняння накладемо графік перехідної функції, побудованої за аналітичною залежністю та графік імпульсної перехідної функції знайденої числовим методом в одній системі координат. Для числового знаходження реакції на використаємо функцію 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. Аналітичні вирази для визначення реакції системи на вхідний синусоїдальний сигнал () за допомогою інтеграла згортки
Реакцію системи (аналітичний її вираз) на вхідний сигнал знайдемо за допомогою iнтеграла згортки:
.
Порівняємо графіки, отримані за допомогою аналітичних виразів та з допомогою функції 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)=(вих-(вх) в усталеному режимі роботи.
Знайдемо амплітудочастотну та фазочастотну характеристику системи, зведеної відносно рівня в другій ємності:
.
Побудуємо амплітудно-частотні та фазо-частотні характеристики системи, використовуючи зовнішню функцію 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
Побудова структурної схеми моделі.
Структурна схема моделі, яка відображає причинно-наслідкові зв’язки між входами та виходами системи:
Рис. 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%
Оголошення від адміністратора