%% prog 002 tubo fuente en un extremo cerrado otro extremo clear all; close all; clc; %este programa calcula los valores propios y vectores propios %los primeros asociados a las frecuencias naturales o de resonancia %y los segundos a la distribución de la presion sonora en el tubo cuando %hay un sonido cuya frecuencia coincide con las frecuencias de resonancia %% inicializacion L = 1.0; %longitud del tubo A = 1e-4; %area de la seccion transversal ro = 1.18; %densidad del aire c = 344.0; %velocidad del sonido Pl = 1.0; %presion sonora en x = -L/2 es decir primer nodo Ne = 100; %numero de elementos N = Ne + 1; %numero de nodos ax = L/(2*Ne); %longitud de elemento %% matriz de masa M = zeros(N,N); %inicializacion de la matriz de masa con ceros %incializacion M(1 ,1 ) = 2; M(N-1,N ) = 1; M(N ,N-1) = 1; M(N , N) = 2; %ciclo principal for n = 2:N-1 M(n ,n) = 4; %diagonal principal M(n ,n-1) = 1; %sub diagonal M(n-1, n) = 1; %super diagonal end M = (A*ax)/(3*c^2)*M; %eliminamos primera fila primera columna M1 = M(2:N,2:N); %nueva matriz de masa FM1 = M(2:N,1); %primera columna de la matriz de masa / fuerza reactiva %% matriz de rigides K = zeros(N,N); %inicializacion de la matriz de rigidez con ceros %incializacion K(1 ,1 ) = 1; K(N-1,N ) = -1; K(N ,N-1) = -1; K(N , N) = 1; %ciclo principal for n = 2:N-1 K(n ,n) = 2; %diagonal principal K(n ,n-1) = -1; %sub diagonal K(n-1, n) = -1; %super diagonal end K = A/(2*ax)*K; %eliminamos primera fila primera columna K1 = K(2:N,2:N); %nueva matriz de rigidez FK1 = K(2:N,1); %primera columna de la matriz de rigidez / fuerza reactiva %% matriz de amortiguamiento zr = ro*c*(1e-2 + j*1e-3); C = zeros(N,N); %inicializacion de la matriz de amortiguamiento con ceros C(N,N) = ro*A/zr; %eliminamos primera fila primera columna C1 = C(2:N,2:N); %nueva matriz de amortiguamiento FC1 = K(2:N,1); %primera columna de la matriz de amortiguamiento / fuerza reactiva % %% problema de valores y vectores propios % % %problema de valores propios % [PHI1,OMG1] = eig(K1,M1); % % %concatenamos una fila de ceros en x = L/2 por la ausencia de fuente % PHI1 = [zeros(1,N-1);PHI1]; % % %frecuencias naturales o de resonancia FEM % wn1 = sqrt(diag(abs(OMG1))); % fnFEM1 = wn1/(2*pi); % % %% frecuencias naturales o de resonancia teoricas TEO % fnTEO1 = zeros(N-1,1); % erro1 = zeros(N-1,1); % nn = 1:N-1; % % for n = 1:N-1 % fnTEO1(n) = ((2*n-1)*c)/(4*L); % erro1(n) = 100*abs(fnFEM1(n) - fnTEO1(n))/fnTEO1(n); % end % % TablaErro1 = [fnFEM1,fnTEO1,erro1]; % % figure(1) % plot(fnFEM1,fnTEO1,'b -o'); % title('frecuencias naturales FEM vs frecuencias naturales TEO'); % xlabel('f_{FEM} (Hz)'); % ylabel('f_{TEO} (Hz)'); % legend('f_n'); % axis([min(fnFEM1),max(fnFEM1),min(fnTEO1),max(fnTEO1)]); % grid on; % box on; % % figure(2) % plot(nn,erro1,'b -o'); % title('error frecuencias naturales FEM vs frecuencias naturales TEO'); % xlabel('n'); % ylabel('error (%)'); % legend('error'); % axis([min(nn),max(nn),0,max(erro1)]); % grid on; % box on; % % %% distribucion de presion sonora en cada frecuencia de resonancia % % dx = L/Ne; %longitud elemento % xcorde = (0:N-1)*dx - L/2; %cordenadas en eje x % x = xcorde; % % figure(3) % subplot(2,2,1) % plot(xcorde, PHI1(:,1)/max( abs( PHI1(:,1) ) ),'b'); % title('primer modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_1'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,2) % plot(xcorde, PHI1(:,2)/max( abs( PHI1(:,2) ) ),'b'); % title('segundo modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_2'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,3) % plot(xcorde, PHI1(:,3)/max( abs( PHI1(:,3) ) ),'b'); % title('tercer modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_3'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,4) % plot(xcorde, PHI1(:,4)/max( abs( PHI1(:,4) ) ),'b'); % title('cuarto modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_4'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % % figure(4) % subplot(2,2,1) % plot(xcorde, PHI1(:,5)/max( abs( PHI1(:,5) ) ),'b'); % title('quinto modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_5'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,2) % plot(xcorde, PHI1(:,6)/max( abs( PHI1(:,6) ) ),'b'); % title('sexto modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_6'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,3) % plot(xcorde, PHI1(:,7)/max( abs( PHI1(:,7) ) ),'b'); % title('septimo modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_7'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,4) % plot(xcorde, PHI1(:,8)/max( abs( PHI1(:,8) ) ),'b'); % title('octavo modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % legend('PHI_8'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; %% resolucion del sistema de ecuaciones en el dominio de la frecuencia %resolveremos el sistema de ecuaciones [-w^2*M1 + j*w*C1 + K1]*P1 = F1 %las frecuencias van desde f = 1 (Hz) hasta 1000 (Hz) %debemos ahora recorda que en el primer nodo la presion sonora P(1) = Pl %coordenadas dx = L/Ne; %longitud elemento xcorde = (0:N-1)*dx - L/2; %cordenadas en eje x x = xcorde; %vector de frecuencias fini = 0; finc = 1; ffin = 1000; f = fini:finc:ffin; w = 2*pi*f; NF = length(f); %prelocalizamos la matriz de solucion donde la columna representa frecuencia %y la fila posicion P1 = zeros(N-1,NF); for nf = 1:NF F1 = w(n)^2*FM1*Pl - j*w(n)*FC1 - FK1*Pl; P1(:,nf) = (-w(nf)^2*M1 + j*w(n)*C1 + K1)\F1; end %colocamos una fila de condicion de contorno primer nodo P(1,1:NF) = Pl = 1 Pini(1,1:NF) = Pl; P1 = [Pini;P1]; %% resultados [Fmesh,Xmesh] = meshgrid(f,x); figure(5) surf(Fmesh,Xmesh,abs(P1),'FaceColor','interp','EdgeColor','none','FaceLighting','phong'); title('modulo de la presion sonora - tubo cerrado - FEM'); xlabel('f (Hz)'); ylabel('x (x)'); zlabel('|P(f,x)| (Pa)'); axis([0,f(NF),-L/2,L/2]); grid on; box on; colormap jet colorbar % %comparamos a una frecuencia fija f = 500 (Hz) % Pf500FEM = P1(:,501); % % %presion teorica % f500 = 500; % k500 = 2*pi*f500/c; % % Pf500TEO = Pl*cos(k500*(x-L/2))/cos(k500*L); % % figure(6) % plot(xcorde,Pf500FEM,'b', xcorde,Pf500TEO,'r -o'); % title('comparación presión sonora f = 500 (Hz) FEM vs TEO'); % xlabel('x (m)'); % ylabel('p(x) (Pa)'); % legend('p(x)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; %presion sonora en función de la frecuencia en x = L/2 PL2 = abs(P1(Ne,:)); LpL2 = 20*log10(PL2/2e-5/sqrt(2)); figure(7) plot(f,LpL2,'b'); title('nivel de presion sonora en x = L/2'); xlabel('f (Hz)'); ylabel('p(f) (Pa)'); legend('p(f)'); axis([fini,ffin,0,1.1*max(LpL2)]); grid on; box on;