% Método de las características % Este programa no considera la aparición de cavitación clc clear all close all % Datos Lgp = 3*154.8 ; % longitud de la galería de presión(m) Ltf = 1197.25-Lgp ; % longitud de la tubería forzada(m) Ltdif = 20 ; % longitud del tubo difusor L= Lgp + Ltf ; landa = 0.02 ; % coeficiente de fricción del conducto D = 4.4 ; % diámetro del conducto (m) % v_0 = 1 ; % m/s c_0 = 1300 ; % m/s N = 50 ; inc_t = 1.0/N ; % Incremento de tiempo tfin = 20.0 ; % tfin = 0.1 ; t_p = (0:1.0/N:tfin); lt=length(t_p); % = tfin*N + 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % nuevos parámetros para imponer el cierre no instantáneo t_c = 1; % tiempo de cierre de la válcula t_iv = L/c_0; %tiempo de ida y vuelta de la onda A_e = (pi*D^2)/4; % Área de entrada de la válvula A_s = A_e ; % Área de salida de la válvula % a_0 = A_s/A_e ; % Constante adimensionales para t=0 % a_0 = 0.01129200682 ; % Para que salga finalmente v_0 próximo a 1 H_g = 400 ; g = 9.81 ; % v_0 = sqrt((2*g*H_g)/(1/(a_0^2)+landa*L/D)); % n = 0.5; % (n<1/2) Cierre muy rápido, siempre aparece compresibilidad % n = 1.0; % Ley de cierre lineal % n = 2.0; % Cierre lento % a = a_0*((t_c-t_p)/t_c).^n ; % Ley de cierre % v_0 = sqrt((2*g*H_g)/(1/(a(1)^2)+landa*L/D)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x_p = (0:1.0/N:1); lx=length(x_p); % = N + 1 rho_0 = 1000 ; % [kg/m3] g = 9.81 ; % [m/s2] H_g = 400; % [m] p_a = 1*10^5 ; % [Pa] alpha = asin(H_g/(Lgp+Ltf+Ltdif)); x = x_p*(Lgp+Ltf+Ltdif); z = (Lgp+Ltf+Ltdif-x)*sin(alpha); landa = 0.02 ; % coeficiente de fricción del conducto D = 4.4 ; % diámetro del conducto (m) Agp = pi/4*D*D ; % área del conducto a la entrada (m^2) Atf = pi/4*D*D ; % área del conducto a la entrada (m^2) % Atdif = 4.5*Atf ; % área del conducto a la salida (m^2) % C1ch = ( landa*Lgp/D + 1.5 )/2/9.81/Agp/Agp ; % CchT = ( landa*Ltf/D + 1.5 )/2/9.81/Atf/Atf ; Ctdif = ( landa*Ltdif/D + 0.5 )/2/9.81/Atf/Atf ; % v_E = v_p*v_0; % H=-103,499 + 22,19·Q - 0,0785307·Q\u2\ B0= -103.499 ; B1= 22.19 ; B2= - 0.0785307 ; % Q = v_p*v_0*Atf ; % HT = B0 + B1*Q + B2*Q.*Q ; % HS = Ctdif*(Atf^2/Atdif^2)*(v_p*v_0).*(v_p*v_0); % HE = HS + HT ; z_E = H_g - (Lgp+Ltf)*sin(alpha); % p_E = (HE - z_E - ((v_p*v_0).*(v_p*v_0))/(2*g))*rho_0*g - p_a; % p_Ep = (p_E+rho_0*g*z_E-p_a-rho_0*g*H_g)/(p_0); t_a = t_c ; t = (0:1.0/N:t_c); % kv = 1 + tanh( (t-t_c*N) / t_a ); C_0=70; o=0.4; % o=1e-6; kv=C_0./(t_c*N-t*N).^o - C_0./(t_c*N).^o; C=0; k = kv; figure(1) clf plot(t,k) % axis([0 200 -1 2*H_g]) xlabel('{\it t}','Fontsize',15) ylabel('{\it k_v}','Fontsize',15) print -deps2 kv_vs_t % Condiciones de iniciales para t=0 (primera fila) y para toda x (todas las columnas) c_0 = 1300 ; % % v_0 = sqrt((2*g*H_g)/(1/(a_0^2)+landa*L/D)); % c_a = (Ctdif+B2+C1ch+CchT)*Atf^2; % c_b = B1*Atf; % % c_c = B0-rho_0*g*H_g; % c_c = B0-H_g; % v_0 = (-c_b + sqrt(c_b^2-4*c_a*c_c))/(2*c_a); % p_0 = rho_0*v_0*c_0; % v_0 = sqrt((2*g*H_g)/(1/(a_0^2)+landa*L/D)); c_a = (Ctdif+B2)*Atf^2+0.5*(-1+landa*L/D); % c_a2 = (Ctdif+B2)*Atf^2+0.5*(1+landa*L/D) c_b = B1*Atf; % c_c = B0-rho_0*g*H_g; c_c = (B0-H_g); v_0 = (-c_b + sqrt(c_b^2-4*c_a*c_c))/(2*c_a) p_0 = rho_0*v_0*c_0; % p_p(1,1:lx) = 0 ; % p_p = rand(1,lx); for i = 1:lx % p_p(1,i) = -(landa*(L+Ltdif)/(2*D))*(v_0/c_0)*i/lx ; p_p(1,i) = -(landa*(L)/(2*D))*(v_0/c_0)*i/lx ; end v_p(1,1:lx) = 1 ; % A continuación, se resuelven el resto de los puntos Q = v_0*Atf ; HT = B0 + B1*Q + B2*Q*Q ; HS = Ctdif*Q*Q; HE = HS + HT ; p_primaE1 = (HE - H_g - v_0^2/2/g)/(v_0*c_0/g) p_primaE2 = p_p(1,lx) A = [1 1 ; 1 -1]; ca = []; cb = []; cc = []; for i = 2:lt j=1; p_p(i,j) = 0 ; v_p(i,j) = v_p(i-1,j+1) - p_p(i-1,j+1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j+1)*abs(v_p(i-1,j+1)) + p_p(i,j) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% j=lx; if i <= t_c*N coefa = ((Ctdif+B2+k(i))*(Atf^2)*(v_0^2)-(v_0^2)/(2*g))*(rho_0*g/p_0); ca = [ca coefa]; coefb = 1+(B1*Atf*v_0)*(rho_0*g/p_0); cb = [cb coefb]; coefc = -v_p(i-1,j-1) - p_p(i-1,j-1) + ((landa*L*v_0*inc_t)/(2*D*c_0))*(v_p(i-1,j-1)^2)+(B0-H_g)*(rho_0*g/p_0); cc = [cc coefc]; % % c_a = (Ctdif+B2)*Atf^2+0.5*(-1+landa*L/D); % c_b = B1*Atf; % c_c = (B0-H_g); v_p(i,j) = (-coefb + sqrt(coefb^2-4*coefa*coefc))/(2*coefa); % Q = v_p(i,j)*v_0*Atf ; % HT = B0 + B1*Q + B2*Q*Q ; % HS = Ctdif*Q*Q; % HE = HS + HT ; % p_E = (HE - z_E - ((v_p(i,j)*v_0)*(v_p(i,j)*v_0))/(2*g))*rho_0*g + p_a; % p_p(i,j) = (p_E+rho_0*g*z_E-p_a-rho_0*g*H_g)/(p_0); % p_p(i,j) = v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)) - v_p(i,j); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% else v_p(i,j) = 0 ; % p_p(i,j) = v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)) - v_p(i,j); end p_p(i,j) = v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)) - v_p(i,j); for j = 2:lx-1 b(1,1) = v_p(i-1,j-1) + p_p(i-1,j-1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j-1)*abs(v_p(i-1,j-1)); b(2,1) = v_p(i-1,j+1) - p_p(i-1,j+1) - landa*L/(2*D)*(v_0/c_0)*inc_t*v_p(i-1,j+1)*abs(v_p(i-1,j+1)); x=A\b; v_p(i,j) = x(1); p_p(i,j) = x(2); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % tcoef = t_p(1:(length(t_p)-2)); tcoef = t_p(1:(t_c*N-1)); figure(2) clf plot(tcoef,ca) xlabel('{\it t}','Fontsize',15) ylabel('{\it ca}','Fontsize',15) figure(3) clf plot(tcoef,cb) xlabel('{\it t}','Fontsize',15) ylabel('{\it cb}','Fontsize',15) figure(4) clf plot(tcoef,cc) xlabel('{\it t}','Fontsize',15) ylabel('{\it cc}','Fontsize',15) % Representación gráfica (animación) % m = size(v_p,1) ; figure(5) clf % hv = plot(x,v(1:n),'EraseMode','background'); hp = plot(x_p,p_p(1,:),'EraseMode','xor'); hold on hv = plot(x_p,v_p(1,:),'r','EraseMode','xor'); axis([0 1 -1.2 1.2]) % axis square % box on set(gca,'xtick',[0 0.2 0.4 0.6 0.8 1]) set(gca,'ytick',[-1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2]) % grid xlabel('{\it x/L}','Fontsize',15) ylabel('{\it v/v}_o , {\it p/p}_o','Fontsize',15) title('Presión y velocidad a lo largo del conducto') for j=2:m pause(0.1); % set(hv,'XData',x,'YData',v( n*(j-1)+1 : n*j )) set(hp,'YData',p_p(j,:)) set(hv,'YData',v_p(j,:)) drawnow end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T = 20; % [ºC] p_vapor = 1000*exp(16.573-3988.842/(273 + T - 39.47)) ; % [Pa] % p_vapor_p = (p_vapor + rho_0*g*z - p_a - rho_0*g*H_g)/p_0; for i = 1:lt for j = 1:lx p(i,j) = p_p(i,j)*p_0 - (rho_0*g*z(j) - p_a - rho_0*g*H_g) ; if (p(i,j) e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.


ESTUDIO NUMÉRICO DE LOS FENÓMENOS TRANSITORIOS APLICADOS A UNA CENTRAL HIDRÃULICA DE BOMBEO PURO

: Moreno Haya, Fernando
: Ingeniería Industrial
Contenido del proyecto:
Directorio raíz  >  Anexos  >  2_Con_compresibilidad  >  Cierre_no_inst_con_fricc_turbina.m