МІНIСТЕРСТВО ОСВIТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ “ЛЬВІВСЬКА ПОЛІТЕХНІКА”
РОЗРАХУНКОВО-ГРАФІЧНА РОБОТА №2
ЗНАХОДЖЕННЯ ДИНАМІЧНИХ ХАРАКТЕРИСТИК
ТЕХНОЛОГІЧНОГО ОБ’ЄКТА ШЛЯХОМ РОЗВ’ЯЗУВАННЯ
СИСТЕМИ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ
Львів - 2008
Для застосування схеми числового інтегрування до системи звичайних диференціальних рівнянь (ЗДР) необхідно її розв'язати відносно перших похідних шуканих функцій. В літературі з числових методів традиційно формули числового інтегрування системи двох ЗДР методом Рунге-Кутта записують у вигляді:
де , - праві частини ЗДР, h - рівномірний крок вибору вузлів інтегрування.
Завдання
Для протічної гідравлічної ємності з коротким та довгим трубопроводами відповідно на вході та виході (рис.1) математична модель (1) має вигляд:
(1)
де , , , , .
Знайти розв’язок системи звичайних диференціальних рівнянь (ЗДР) - H(t), Q1(t) (зміну рівня в ємності та витрати рідини в трубопроводі) на інтервалі t([0; 4], c з кроком Δt=0,05 c, якщо
d=0,25 м
r1=0,06 м
L1=90 м
r2=0,06 м
L2=2 м
P1=1500 Па
P2=300 Па
(=0,9
ρ=1e3 кг/м3
g=9,81 м/с2
ν=1e-5 Па
Початкові умови: значення рівня H(0)=0,1 м та витрати Q1(0)=0,001 м3/с.
Розв’язок подати у вигляді графіка.
Розв’язання
Застосуємо цю схему числового інтегрування системи двох ЗДР для системи (1) диференціальних рівнянь. Спочатку розв’яжемо її відносно перших похідних шуканих функцій, тобто
(2)
Позначимо
;
або
.
Тоді схема числового інтегрування системи двох ЗДР (2) матиме вигляд: (3)
Програма числового розв'язування системи двох ЗДР мовою С.
#include <stdio.h>
#include <conio.h>
#include <math.h>
float f1(float t, float H, float Q1)
{float d=0.25,r2=0.06,L2=2,g=9.81,ro=1e3,nu=1e-5;
float S=M_PI*d*d/4,P2=300;
float k2=M_PI*pow(r2,4)/8*nu*L2;
return 1/S*(Q1-(k2*ro*g*H-P2)/ro);
}
float f2(float t, float H, float Q1)
{float r1=0.06,L1=90,dz=0.9,ro=1e3;
float k1=2*M_PI*r1*r1*sqrt(r1/L1*dz), P1=1500;
float A=4*M_PI*pow(r1,3)/dz;
return 1/A*(k1*k1*P1/ro-Q1*Q1);
}
main()
{clrscr();
float t, t0=0, tf=4, n=40;
float h,k1,k2,k3,k4,l1,l2,l3,l4,Hi=0.1,Q1i=0.001;
h=(tf-t0)/n;
for(t=t0;t<=tf;t+=h)
{printf("\n t=%f Hi=%f Q1i=%f",t,Hi,Q1i);
k1=h*f1(t, Hi,Q1i);
l1=h*f2(t, Hi,Q1i);
k2=h*f1(t+h/2, Hi+k1/2, Q1i+l1/2);
l2=h*f2(t+h/2, Hi+k1/2, Q1i+l1/2);
k3=h*f1(t+h/2, Hi+k2/2, Q1i+l2/2);
l3=h*f2(t+h/2, Hi+k2/2, Q1i+l2/2);
k4=h*f1(t+h, Hi+k3, Q1i+l3);
l4=h*f2(t+h, Hi+k3, Q1i+l3);
Hi=Hi+1.0/6*(k1+2*k2+2*k3+k4);
Q1i=Q1i+1.0/6*(l1+2*l2+2*l3+l4);
}
getchar();
return 0;
}
Запишемо в таблицю 1 кожне десяте значення функцій для рівня H(t), м та витрати Q1(t), м3/с.
Таблиця 1
t,c
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
H(t),м
0.1000
0.7132
1.3263
1.9394
2.5525
3.1656
3.7786
4.3916
5.0045
5.6175
Q1(t), м3/с
0.0010
0.0010
0.0010
0.0009
0.0009
0.0009
0.0009
0.0009
0.0009
0.0009
t,c
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
H(t),м
6.2304
6.8433
7.4562
8.0691
8.6820
9.2948
9.9076
10.5205
11.1333
11.7461
Q1(t), м3/с
0.0009
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
t,c
2
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
H(t),м
12.3588
12.9716
13.5843
14.1971
14.8098
15.4226
16.0353
16.6480
17.2607
17.8734
Q1(t), м3/с
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0008
0.0007
t,c
3
3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8
3.9
4
H(t),м
18.4861
19.0987
19.7114
20.3241
20.9367
21.5494
22.1620
22.7747
23.3873
23.9999
24.6126
Q1(t), м3/с
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
Розв’яжемо систему двох ЗДР засобами MatLab.
Створимо файл dataf.m, в якому запишемо вхідні дані та обчислимо коефіцієнти.
%data5
ro=1e3; mu=1e-5; g=9.81; dz=0.9;
d=0.25;
S=pi*d^2/4;
P1=15e3; P2=3e3;
L1=90; r1=0.06;
L2=2; r2=0.06;
A=4*pi*r1^3/dz;
k1=2*pi*r1^2*sqrt(r1/L1/dz);
k2=pi*r2^4/8/mu/L2;
Файл-функція f.m, де записані праві частини системи двох ЗДР, розв’язаних відносно перших похідних, має вигляд:
function y=f5(t,x);
Q1=x(1); H=x(2);
data5;
y=[1/A*(k1^2/ro*P1-Q1^2);
1/S*(Q1-k2*(ro*g*H-P2)/ro)];
Дані та обчислені значення коефіцієнтів передають в файл-функцію f.m шляхом включення файлу даних dataf.m.
Файл main_Matlab.m, в якому задають інтервал часу [t0, tf], на якому інтегрують систему ЗДР, початкові умови (x0=[1.0 0.01]) та записують функцію MatLab ode45(), призначену для розв’язування ЗДР чи системи ЗДР числовим методом має вигляд:
%main5
t0=0;tf=4;x0=[0.001 0.1];
[t,y]=ode45('f5',[t0 tf],x0);
figure(1); subplot(2,1,1);plot(t,y(:,1),'r');grid;xlabel('t,c');ylabel('Q_1,m^3/c');
subplot(2,1,2);plot(t,y(:,2),'r');grid;xlabel('t,c');ylabel('H,m');
В цьому файлі дві останні стрічки призначені для побудови графіків розв’язків системи заданої ЗДР.
Рис.2 Графіки перехідних процесів:
а) зміни рівня в ємності; б) витрати рідини в трубопроводі.
Накладемо графіки, отримані в засобами MatLab та в Сі. Для цього модифікуємо файл main5.m, додавши до нього дані таблиці 1.
%main5
t0=0;tf=4;x0=[1.0 0.01];
[t,y]=ode45('f5',[t0 tf],x0)
ti=[0:0.05:4];
Q1i=[ 0.1000 0.7132 1.3263 1.9394 2.5525 3.1656 3.7786 4.3916 5.0045 5.6175 6.2304 6.8433 7.4562 8.0691 8.6820 9.2948 9.9076 10.5205 11.1333 11.7461 12.3588 12.9716 13.5843 14.1971 14.8098 15.4226 16.0353 16.6480 17.2607 17.8734 18.4861 19.0987 19.7114 20.3241 20.9367 21.5494 22.1620 22.7747 23.3873 23.9999 24.6126];
Hi=[0.1000 0.7132 1.3263 1.9394 2.5525 3.1656 3.7786 4.3916 5.0045 5.6175 6.2304 6.8433 7.4562 8.0691 8.6820 9.2948 9.9076 10.5205 11.1333 11.7461 12.3588 12.9716 13.5843 14.1971 14.8098 15.4226 16.0353 16.6480 17.2607 17.8734 18.4861 19.0987 19.7114 20.3241 20.9367 21.5494 22.1620 22.7747 23.3873 23.9999 24.6126];
figure(1); subplot(2,1,2);plot(t,y(:,2),'k',ti,Q1i,'b*');grid;xlabel('t,c'); ylabel('Q_1,m^3/c');
subplot(2,1,1);plot(t,y(:,1),'k',ti,Hi,'b*');grid;xlabel('t,c'); ylabel('H,m');
Рис.3 Порівняння графіків перехідних процесів, отриманих в MatLab (-) та в Сі (*) для:
а) зміни рівня в ємності; б) витрати рідини в трубопроводі.
.
Висновок. Розв’язки системи ЗДР, отримані засобами MatLab та в Сі збігаються, а значить можна зробити висновок, що знайдені правильно.