method of characteristics, boundary conditions

Question

Kirill T 2016-4-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/278211-method-of-characteristics-boundary-conditions

回答： badria alrashidi 2020-11-22

Good day! Tell who knows. Compiled the code to calculate the wave processes in the hydraulic lines (solve the problem of the method of characteristics). If we consider one cycle (without operator "while"), it's done. When add the operator "while" produces an empty multiplicity. The program code below. The result is that some should be able to see Photo

I tried to realize the next code in matlab, see adede files

if true
  % code
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;
end