Побудова математичної моделі та числове її дослідження

 

 

Конструктивні параметри: L1=100м, r1=0.05м, kB=0.005м , d=1м.

Вхідні величини: P1=2000Па, P2=9500Па,T1=293K, T2=350K.

Керуючі величини: l2=0.4, l3=0.7.

Значення збурення:  .

Параметри стану системи: h,T.

Згідно з рівнянням збереження маси речовини та ввівши деякі припущення (масообмін на границі розділу фаз рідина-повітря відсутній, =const), запишемо диференційне рівняння, що описує зніму рівня в ємності:

 

 ; (1.1)

 

де  - площа ємності ,м ;

h- рівень рідини в ємності, м;

Q- об*ємна витрата, м /с. 

  визначається наступною залежністю:

- для короткого турбулентного трубопроводу з регулюючим органом (P3=0): 

 

Об*ємна витрата Q1 визначається за формулою:

 -для довгого ламінарного трубопроводу:

 

  (1.2)

 

Об*ємна витрата Q2 визначається наступною залежністю:

- для короткого турбулентного трубопроводу з регулюючим органом : 

 

 Для отримання виразу зміни температури рідини в ємності, скористаємося законом збереження тепла. Запишемо:

 

  (1.3)

 

або 

 

  (1.3*)

 

Враховуючи рівняння (1.1) і (1.3*) зробимо заміну   і отримуємо:

  

   (1.4)

 

Об*єднавши рівняння (1.1), (1.2), (1.4) отримаємо наступну систему звичайних диференціальних рівнянь:

 

  (1.5)

 

Ця система звичайних диференційних рівнянь разом із початковими умовами   є математичною моделлю даного об*єкту. 

 

 

1. Знайдемо початкові умови.

Знайдемо розв*язок системи нелінійних рівнянь застосовуючи функцію fsolve . Для цього створимо файл даних yos_kr1.m :

 

%yos_kr1.m

function y=yos_kr1(x);

h=x(1); Q1=x(2); T=x(3);

%-------------------------------

P1=2000;P2=9500; T1=293; T2=350;

L1=100; r1=0.05; kv=0.005; 

l2=0.4; l3=0.7; d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

%-------------------------------

S=pi*d^2/4; 

kl=pi*r1^4/8/L1/nu;

T=L1*kl/pi/r1^2;

Q2=kv*l2*sqrt((P2-ro*g*x(1))/ro);

Q3=kv*l3*sqrt(g*x(1));

y=[(x(2)+kv*l2*sqrt((P2-ro*g*x(1))/ro)-kv*l3*sqrt(g*x(1)))/S;

 (kl*P1/ro-x(2))/T;

 (x(2)*(T1-x(3))+kv*l2*sqrt((P2-ro*g*x(1))/ro)*(T2-x(3)))/S/x(1)];

 

Потім створюємо файл yos_kr11.m:

 

x0=[0.1 0.003 320];

y=fsolve(*yos_kr1*,x0)

 

Результатами виконання файлу є такі значення параметрів стану: 

 

h0=0.617917411421; Q10=0.004908738521; T0=317.530396154693;

 

Знайдемо тепер номінальні значення параметрів стану об*єкту числовим методом, використовуючи функції MatLab. Для того, щоб розв*язати систему нелінійних диференціальних рівнянь (1.5) потрібно створити два файли, в одному з яких (файлі-функції yos_kr2.m) будуть записані праві частини системи диференціальних рівнянь, розв*язані відносно перших похідних, а у другому файлі (yos_kr21.m) буде записана функція MatLab ode45, призначена для розв*язування системи диференціальних рівнянь.

 

Файл-функції yos_kr2.m:

function y=yos_kr2(t,x);

h=x(1); Q1=x(2); T=x(3);

%----------------------------

P1=2400;P2=9500; T1=293; T2=350;

L1=100; r1=0.05; kv=0.005; 

l2=0.4; l3=0.7; d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

%----------------------------

S=pi*d^2/4; 

kl=pi*r1^4/8/L1/nu;

T=L1*kl/pi/r1^2;

Q2=kv*l2*sqrt((P2-ro*g*x(1))/ro);

Q3=kv*l3*sqrt(g*x(1));

%--------------------------------------------------------

y=[(x(2)+kv*l2*sqrt((P2-ro*g*x(1))/ro)-kv*l3*sqrt(g*x(1)))/S;

(kl*P1/ro-x(2))/T;

(x(2)*(T1-x(3))+kv*l2*sqrt((P2-ro*g*x(1))/ro)*(T2-x(3)))/S/h];

 

Файл yos_kr21.m:

x0=[0.617917411421 0.004908738521 317.530396154693];

T=[0 500];

tol=odeset(*abstol*,[1e-20 1e-20 1e-20]);

[t,y]=ode45(*yos_kr2*,T,x0,tol) ;

plot(t,y(:,1));grid; 

ylabel(*h,m*); 

xlabel(*t,sec*);pause;

plot(t,y(:,2));grid;

ylabel(*Q1,mkub/sec*);

xlabel(*t,sec*);pause;

plot(t,y(:,3));grid;

ylabel(*T,K*);

xlabel(*t,sec*);

Результатом роботи цих програм є графіки перехідних процесів параметрів стану досліджуваного об*єкту.

 

Рисунок 1.1 Перехідний процес h(t) з номінального початкового стану рівноваги до нового стану рівноваги, при нанесенні на об*єкт збурення  .

 

Рисунок 1.2 Перехідний процес Q1(t) з номінального початкового стану рівноваги до нового стану рівноваги, при нанесенні на об*єкт збурення  .

 

Рисунок 1.3 Перехідний процес T(t) з номінального початкового стану рівноваги до нового стану рівноваги, при нанесенні на об*єкт збурення  .

 

 

 

Запишемо систему (1.5) у вигляді

 

  (2.1)

де  

 

Лінеаризована система диференційних рівнянь матиме вигляд:

 

  (2.2)

 

або

 

  (2.3)

 

вихідні величини : y1(t)=h (t)

 

y2(t)=T(t)

 

Часткові похідні правих частин системи нелінійних диференційних рівнянь (2.1) беремо по всіх параметрах стану, а також по тих вхідних величинах, відхилення яких від номінального режиму задане в завданні на курсову роботу: 

 

Коефіцієнти   матриці стану системи та вектора   вхідних величин обчислюються наступним чином:

 

В матричній формі лінеаризована система диференціальних рівнянь матиме вигляд: 

 

 , (2.4)

де   - матриця стану системи;

  - вектор вхідних величин;

 ;  .

  - вектор параметрів стану системи;

  - сигнал збурення.

 

 Для обчислення коефіцієнтів власної у матриці системи та вектора вхідних величин складемо програму yos_kr3.m:

 

h0=0.617917411421; Q10=0.004908738521; T0=317.530396154693;

%--------------------------------------------------------

P1=2000;P2=9500; T1=293; T2=350;

L1=100; r1=0.05; kv=0.005; 

l2=0.4; l3=0.7; d=1; nu=0.00001;

dz=0.9; ro=1e3; g=9.81;

S=pi*d^2/4; 

kl=pi*r1^4/8/L1/nu;

T=L1*kl/pi/r1^2;

%----------------------------------

a11=-(kv*l2*g/2/sqrt((P2-ro*g*h0)/ro)+kv*l3*g/2/sqrt(g*h0))/S;

a12=1/S;