Реализация метода характеристик (расчет волновых процессов в гидролиниях)

clear all; clc;
diam = 50*10^-3;%input('Введите внутренний диаметр трубы, мм: ')/1000;
Area = (pi/4)*diam*diam;   %* [m2] *
q0 = 3/(1000*60);%input('Введите значение расхода, л/мин: ')/(1000*60);
v0 = q0/(Area);        % средняя скорость * [m/s] *
len = 800; %input('Введите длину трубы, м: ');
cel = 1000; %input('Введите скорость звуковой волны, м/с: ');
n = 1000; % Число узлов (number of computational sections)
ds = len/n; % шаг по длине (шаг дискритизации пространства - spatial discretization step)
dt = ds/cel; % шаг по времени
h0 = 100; % напор в баке, м
Roughness = 0.0125;
g = 9.806; %Ускорение свободного падения, м/с2
Lambda = 124.528/(power((1/Roughness),2)*power(diam,(1/3)));
tc = 0; %время включения клапана на выходе, с
m = 1; % Показатель степени закона закрытия задвижки
tmax = 15; %Время вычисления
 
 
Vp = 0;
vp = 0;
Hp = 0;
hp = 0;
NegligibleHeadlosses = 0;
     t = 0;
     while t < tmax 
          t = t + dt;
            if NegligibleHeadlosses == 1
                for i = 1:n+1
                    h(i) = h0;
                    v(i) = v0;
                    Lambda = 0;
                end;
            else
                for i = 1:n+1
                    hf = Lambda*ds*v0*v0/(2*g*diam);
                    h(i) = h0-i*hf;
                    v(i) = v0;
                end;
            end;
        for i = 2:n
          vp(i) = 0.5*(v(i-1)+v(i+1)+(g/cel)*(h(i-1)-h(i+1))-(Lambda*dt/(2*diam))*(v(i-1)*abs(v(i-1))+v(i+1)*abs(v(i+1))));
          hp(i) = 0.5*(h(i-1)+h(i+1)+(v(i-1)-v(i+1))/(g/cel)-(Lambda*dt/(2*diam))*(v(i-1)*abs(v(i-1))-v(i+1)*abs(v(i+1)))/(g/cel));
        end;
 %(* Reservoir positioned upstream *)
     hp(1) = h0;
     vp(1) = v(2)+g*(h(1)-hp(2))/cel-Lambda*v(2)*abs(v(2))*dt/(2*diam);     
 
     
%(* Valve positioned downstream *)
     if t < tc
          tau = power((1-t/tc),m);
          vp(n+1) = v0*tau;
          hp(n+1) = h(n)-cel*(vp(n+1)-v(n))/g-Lambda*cel*v(n)*abs(v(n))*dt/(2*g*diam);
     else
          tau = 0;
          vp(n+1) = 0;
          hp(n+1) = h(n)+cel*v(n)/g-Lambda*cel*v(n)*abs(v(n))*dt/(2*g*diam);
     end;
     
          for i = 1:n+1
%                v(end+1,: ) = vp(i);
%                h(end+1,: ) = hp(i);
          end;
%           Vp = vp(i) + Vp;
%           Hp = h(i) + Hp;
     end;
   
t = 1:dt:tmax;     
figure;
plot(t,h,'r');
axis([0 tmax 0 220]);
grid on;
hold on;

if NegligibleHeadlosses == 1
k=1:n+1;
h= h0.*k;
v = v0.*k;
Lambda = 0;
else 
k=1:n+1;
hf = (Lambda*ds*v0*v0/(2*g*diam)).*k;
h = h0-k.*hf;
v = v0.* ones(1,n+1);
end;

figure
plot(x,y)
hold on
plot(x,y)