%% Prog010FEMTuboFuenteCerradoImpRadConMontaje clearvars; close all; clc; %este es un programa de elementos finitos que representa un tubo abierto en un extremo %con fuete y cerrado en el otro extremo %en ambos extremos calcularemos las frecuencias naturales (valores propios) %y la distribucion de la presion sonora (MNV) dentro del tubo a dichas %frecuencias (vectores propios) %para despues determinar la respuesta forzada %% datos inciales Ne = 100; %numero de elementos NnD = 0; %numero de DOF s con condiciones drichelet N = Ne+1; %numero de DOFs - funciones de interpolacion Nt = Ne+1-NnD; %numero de DOFs - funciones de interpolacion con drichelet Nnod = 2; %numero de nodos elemento Ncord = 2; %numero de coordenadas por elemento Ndof = 1; %numero de dof por nodo L = 1.0; %longitud tubo m ax = L/(2*Ne); %logitud de 1/2 elemento a = 0.01; %lado/radio del tubo A = pi*a^2; %area de la seccion transversal del tubo A c = 344; %velocidad del sonido ro = 1.18; %densidad del aire B = c^2*ro; %rigidez de compresibilidad adiabática de volumen P0 = 2e-5; %recordemos que este es un modelo de ondas planas %esto se cumple cuando la longitud de onda es mucho mayor que las %dimensiones que forman el area de la seccion transversal %se utilizara matrices elementales %% matriz de cordenadas %inicio Xcord = zeros(N,3); %coordenadas deltax = L/Ne; xcorde = transpose( (0:N-1)*deltax - (L/2) ); Xcord(:,1) = xcorde; %% matriz de conectividad %inicio Xconec = zeros(Ne,Ncord); for n = 1:N-1 Xconec(n,1) = n; Xconec(n,2) = n+1; end %% matriz de ID %inicio XID = zeros(N,Ndof); %nodos con condicion de drichelet (ninguno en este caso) nodDr = ones(N,Ndof); 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 - rigidez elemental y global %inicializacion matriz de masa global M = zeros(Nt,Nt); %matriz de masa elemental Me = [2 1;1 2]; Me = (A*ax)/(3*c^2)*Me; %inicializacion matriz de rigidez global K = zeros(Nt,Nt); %matriz de 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; %% correccion de fuente M1 = M(2:Nt,2:Nt); K1 = K(2:Nt,2:Nt); C1 = C(2:Nt,2:Nt); % %% valores y vectores propios % % %la función eig() determina los valores y vectores propios de la matriz % [PHI1,LAM1] = eig(K1,M1); % % %incorporamos la condicion de contorno de tubo abierto con la fuente % %apagada para ello incorporamos una fila llena de ceros % PHI1 = [zeros(1,N-1);PHI1]; % % %frecuencias naturales elementos finitos FEM % f_FEM1 = diag( real(sqrt(LAM1)))/(2*pi); % nn = 1:length(f_FEM1); % % %frecuencias naturales teóricas TEO % f_TEO1 = zeros(N-1,1); % for n = 1:N-1 % f_TEO1(n) = (2*n-1)*c/(4*L); % end; % % %error % err1 = 100*abs(f_FEM1 - f_TEO1)./f_TEO1; % % %tabla % f_FEM1_f_TEO1_err1 = [f_FEM1, f_TEO1, err1]; % % %% figuras frecuencias valores propios % % figure(1) % plot(f_TEO1,f_FEM1,'b -o'); % title('frecuencias f_{TEO} vs f_{FEM}'); % xlabel('f_TEO (Hz)'); % ylabel('f_TEO (Hz)'); % axis([min(f_TEO1),max(f_TEO1),min(f_FEM1),max(f_FEM1)]); % grid on; % box on; % % figure(2) % plot(nn,err1,'b -o'); % title('porcentaje de error frecuencias f_{TEO} vs f_{FEM}'); % xlabel('n'); % ylabel('error (%)'); % grid on; % box on; %% figuras distribución de presion sonora modal %coordenadas nodales x = zeros(N,1); for n = 1:N x(n) = -(L/2) + (n-1)*L/Ne; end % figure(3) % % subplot(2,2,1) % plot(x, PHI1(:,1)/max( abs(PHI1(:,1) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(1))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,2) % plot(x, PHI1(:,2)/max( abs(PHI1(:,2) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(2))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,3) % plot(x, PHI1(:,3)/max( abs(PHI1(:,3) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(3))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,4) % plot(x, PHI1(:,4)/max( abs(PHI1(:,4) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(4))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % figure(4) % % subplot(2,2,1) % plot(x, PHI1(:,5)/max( abs(PHI1(:,5) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(5))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,2) % plot(x, PHI1(:,6)/max( abs(PHI1(:,6) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(6))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,3) % plot(x, PHI1(:,7)/max( abs(PHI1(:,7) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(7))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; % % subplot(2,2,4) % plot(x, PHI1(:,8)/max( abs(PHI1(:,8) ) ) ,'b') ; % title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(8))] ); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2,L/2,-1.5,1.5]); % grid on; % box on; %% solucion forzada dominio de la frecuencia metodo directo %resolveremos el sistema de ecuaciones [-w^2*M1 +j*w*C1 + K1]*P1(w) = F1(w) %donde w = 2*pi*f y f = 1,2,3,...,1000 (Hz) %la pesion sonora en x = -L/2 es PL = exp(-j*w*t) con PL = 1 (Pa) %componentes del vector de fuerzas PL = 1e-3; %amplitud de presion sonora Pm = M(2:N,1); %primera columna de la matriz de masa Pc = C(2:N,1); %primera columna de la matriz de masa Pk = K(2:N,1); %primera columna de la matriz de rigidez %frecuencias fmin = 1; fmax = 1000; finc = 1; f = fmin:finc:fmax; %vector de frecuencias w = 2*pi*f; %vector de frecuencias angulares NF = length(f); %prelocalizamos la solucion como matriz %la primera columna tiene la solución para f = 1 (Hz) %la segunda columna tiene la solución para f = 2 (Hz) %................................................... %la ultima columna tiene la solución para f = 1000 (Hz) 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); %impedancia arbitraria %ZL = ro*c*( 0.1 + j*0.01 ); 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 k = 2*pi*f(500)/c; p = 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(1) surf(Fmesh,Xmesh,abs(P1),'FaceColor','interp',... 'EdgeColor','none',... 'FaceLighting','phong') title('Modulo Presion Sonora - Tubo Impedancia Radiacion - 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(2) surf(Fmesh,Xmesh,20*log10((abs(P1) + P0)/sqrt(2)/2e-5),'FaceColor','interp',... 'EdgeColor','none',... 'FaceLighting','phong') title('Nivel de Presion Sonora - Tubo Impedancia Radiacion - 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) %p(x,t) = PL*sin[k*(x-L/2)]/sin(k*L) figure(3) plot(x,abs(p),'b',x,abs(P1(:,500)),'r -o'); title(['comparacion solucion TEO vs FEM frecuencia f = ',num2str(f(500))] ); xlabel('x (m)'); ylabel('p (Pa)'); legend('p_{TEO}(x,f)','p_{FEM_DIR}(x,f)'); axis([-L/2,L/2,0.0,1.5*max(p)]); grid on; box on; figure(4) plot(f,20*log10(P1(1,:)/sqrt(2)/P0),'b',f,20*log10(P1(Nt,:)/sqrt(2)/P0),'r'); title('niveles de presion sonora al incio y al final del tubo'); xlabel('f (Hz)'); ylabel('L_p (dB)'); legend('L_{p1}(f)','L_{p2}(f)'); grid on; box on; figure(5) plot(f,20*log10( P1(1,:)./P1(Nt,:)),'b'); title('reduccion de ruido NR'); xlabel('f (Hz)'); ylabel('NR(f) (dB)'); legend('NR(f)'); grid on; box on; figure(6) plot(f,20*log10( P1(Nt,:)./P1(1,:)),'b'); title('respuesta de frecuencia'); xlabel('f (Hz)'); ylabel('H(f) (dB)'); legend('H (f)'); grid on; box on;