%% prog 010 tubo cerrado fuente en un extremoimpedancia de radiacion otro con montaje clearvars; 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 = 0.01; %lado/radio del tubo A = pi*a^2; %area de la seccion transversal del tubo A ro = 1.18; %densidad del aire c = 344.0; %velocidad del sonido Ne = 100 ; %numero de elementos N = Ne + 1; %numero de nodos Ncord = 2; %numero de cordenadas por elemento Nnod = 2; %numero de nodos por elemento Ndof = 1; %numero de grados de libertad por nodo ax = L/(2*Ne); %longitud de elemento PL = 1.0; %fuente de presion en x = -L/2 %% matriz de coordenadas %inicializacion Xcord = zeros(N,3); %coordenadas en x deltaX = L/Ne; xcorde = transpose( (0:N-1)*deltaX - L/2 ); Xcord(:,1) = xcorde; %% matriz de conectividad Xconec = zeros(Ne,Ncord); for n = 1:Ne Xconec(n,1) = n; Xconec(n,2) = n + 1; end; %% matriz de condiciones de Dirichelet XID = zeros(N,Ndof); %matriz de condiciones de Dirichelet nodDr = ones(N,Ndof); %nodo final y quinto nodo abierto salen cosas raras % nodDr(N,1) = 0; % nodDr(5,1) = 0; cont = 1; %contador for n = 1:N if nodDr(n) == 0 XID(n,1) = 0; else XID(n,1) = cont; cont = cont + 1; end; end; %% matriz de masa y rigidez M = zeros(N,N); %inicializacion de la matriz de masa global con ceros K = zeros(N,N); %inicializacion de la matriz de rigidez global con ceros %matriz masa elemental Me = [ 2, 1 ; 1, 2]; Me = (A*ax)/(3*c^2)*Me; %matriz rigidez elemental Ke = [ 1,-1 ; -1, 1]; Ke = A/(2*ax)*Ke; %loop principal montaje for nel = 1:Ne for noi = 1:Nnod for ngi = 1:Ndof for noj = 1:Nnod for ngj = 1:Ndof NN = XID( Xconec(nel,noi) , ngi ); MM = XID( Xconec(nel,noj) , ngj ); if (NN ~= 0) && (MM ~= 0) fila = Ndof*(Nnod-1)*(noi - 1) + ngi; colu = Ndof*(Nnod-1)*(noj - 1) + ngj; M(NN,MM) = M(NN,MM) + Me( fila, colu ); K(NN,MM) = K(NN,MM) + Ke( fila, colu ); end; end; end; end; end; end; %% matriz de amortiguamiento C = zeros(N,N); C(N,N) = 1; %matriz de amortiguamiento con condicion de contorno de drichelet no nula C1 = C(2:N,2:N); %% fuente sonora %en este caso implica eliminar primera fila y primera columna del sistema %de ecuaciones %usar parte de la primera fila / columna de cada matriz para armar la %fuerza reactiva %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 %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 %% 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 % % erro(1) = 0; % 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 dx = L/Ne; %longitud elemento xcorde = (0:N-1)*dx - L/2; %ccordenadas en eje x % 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 + 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 %vector de frecuencias fini = 0; finc = 1; ffin = 1000; f = fini:finc:ffin; w = 2*pi*f; NF = length(f); %componentes del vector de fuerzas PL = 1; %amplitud de presion sonora Pm = M(2:N,1); %primera columna de la matriz de masa Pk = K(2:N,1); %primera columna de la matriz de rigidez Pc = C(2:N,1); %primera columna de la matriz de amortiguamiento P1 = zeros(N-1,NF); %solucion for n = 1:NF %impedancia de radiacion %ZL = ro*c*( 1-2*besselj(1,2*k*a) + 2*j*struve(1,2*k*a) ) k = w(n)/c; ka2 = 2*k*a; RL = 1 - 2*besselj(1,ka2)/ka2; XL = 2*struve(1,ka2)/ka2; ZL = ro*c*( RL + j*XL); C1(N-1,N-1) = ro*A/ZL; %vector de fuerza dependiente de la frecuencia F1 = w(n)^2*PL*Pm - j*w(n)*Pc - PL*Pk; %usamos A\b que es mejor que inv(A)*b P1(:,n) = (-w(n)^2*M1 +j*w(n)*C1 + K1)\F1; end; %incorporamos la solucion PL para x = -L/2 P1 = [PL*ones(1,NF);P1]; %% comparación resultados teóricos y elementos finitos x = xcorde; k = 2*pi*f(600)/c; pcer = PL*cos(k*(L/2 - x))/cos(k*L); pabr = PL*sin(k*(L/2 - x))/sin(k*L); %grafico espacio - frecuencia - presion sonora %distribucion del valor absoluto de la presion sonora ontenida usando FEM %en funcion de la posicion "x" y de la frecuencia "f" [Fmesh,Xmesh] = meshgrid(f,x); figure(5) surf(Fmesh,Xmesh,abs(P1),'FaceColor','interp',... 'EdgeColor','none',... 'FaceLighting','phong') title('Modulo Presion Sonora - Tubo Cerrado - FEM'); xlabel('f (Hz)'); ylabel('x (m)'); zlabel('|P| (Pa)') axis([0,f(NF),-L/2,L/2]); grid on; box on; colormap jet colorbar; figure(6) surf(Fmesh,Xmesh,20*log10(abs(P1)/sqrt(2)/2e-5),'FaceColor','interp',... 'EdgeColor','none',... 'FaceLighting','phong') title('Nivel de Presion Sonora - Tubo Cerrado - FEM'); xlabel('f (Hz)'); ylabel('x (m)'); zlabel('|P| (Pa)') axis([0,f(NF),-L/2,L/2]); grid on; box on; colormap jet colorbar; %compararemos con la solucion exacta para la frecuencia de f = 1000 (Hz) %pcer(x,t) = PL*cos[k*(L/2-x)]/cos(k*L) %pabr(x,t) = PL*sin[k*(L/2-x)]/sin(k*L) figure(7) plot(x,abs(pcer),'b',x,abs(pabr),'k',x,abs(P1(:,600)),'r -o'); title(['comparacion solución TEO vs FEM frecuancia f = ',num2str(f(600))] ); xlabel('x (m)'); ylabel('p (Pa)'); legend('p_{TEO_{CER}}(x,f)','p_{TEO_{ABR}}(x,f)','p_{FEM_{DIR}}(x,f)'); axis([-L/2,L/2,-0.1*max(abs(pcer)),2.0*max(abs(pcer))]); grid on; box on;