МІН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 та в Сі збігаються, а значить можна зробити висновок, що знайдені правильно.