method of characteristics, boundary conditions

16 views (last 30 days)
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

Answers (1)

badria alrashidi
badria alrashidi on 22 Nov 2020
hi I need code to apply in matlab pdes method of Characteristics

Categories

Find more on Function Creation in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!