%% prog 005 FEM tubo fuente en extremo 1 cerrado 2 impedancia de radiacion sin montaje clearvars; close all; clc; %este programa no calcula de forma directalos valores propios %(frecuencias naturales) y vectores propios (formas de vibracion) de un %tubo cerrado en ambos extremos de seccion transversal constante A %longitud L a una velocidad del sonido constante c %inicialmente el numero de elementos es 4 pero será modificado a un numero %variable de elementos %% inicializacion Ne = 100; %numero de elementos N = Ne + 1; %numero de nodos / grados de libertad / DOF L = 1.0; %longitud del tubo en metros (m) ax = L/(2*Ne); %mitad de la longitud de un elemento A = 0.0001; %area de la sección transversal (m^2)( 1.0e-4) = 1 X 10^(-4) ) c = 343; %velocidad del sonido ro = 1.18; %densidad del aire rd = sqrt(A/pi); %radio de tubo; %% matriz de masa %inicializacion M = zeros(N,N); %matriz incializada con ceros %en una primera etapa colocaremos de forma directa todos los valores en la %matriz %M = [2 1 0 0 0; % 1 4 1 0 0; % 0 1 4 1 0; % 0 0 1 4 1; % 0 0 0 1 2]; %utilizaremos la estructura de la matriz para generalizar resultados para %tubos de seccion cosntante donde todos los elementos tienen el mismo %tamaño %proceso incial M(1 , 1 ) = 2; M(N-1, N ) = 1; M(N , N-1) = 1; M(N , N ) = 2; %iteracion principal desde la 2da fila/columna a la penultima fila/columna for n1 = 2:N-1 M(n1 ,n1 ) = 4; %diagonal llena de 4 M(n1-1,n1 ) = 1; %super diagonal llena de 1 M(n1 ,n1 - 1) = 1; %sub diagonal llena de 1 end; %factores de ajuste M = (A*ax)/(3*c^2)*M; %% matriz de rigidez %inicializacion K = zeros(N,N); %en una primera etapa colocaremos de forma directa todos los valores en la %matriz %K = [ 1 -1 0 0 0; % -1 2 -1 0 0; % 0 -1 2 -1 0; % 0 0 -1 2 -1; % 0 0 0 -1 1]; %utilizaremos la estructura de la matriz para generalizar resultados para %tubos de seccion cosntante donde todos los elementos tienen el mismo %tamaño %proceso incial K(1 , 1 ) = 1; K(N-1, N ) = -1; K(N , N-1) = -1; K(N , N ) = 1; %iteracion principal desde la 2da fila/columna a la penultima fila/columna for n1 = 2:N-1 K(n1 ,n1 ) = 2; %diagonal llena de 2 K(n1-1,n1 ) = -1; %super diagonal llena de -1 K(n1 ,n1 - 1) = -1; %sub diagonal llena de -1 end; %factores de ajuste K = A/(2*ax)*K; %% matriz de amortiguamiento %impedancia de terminacion tubo cerrado %Z = 10e3 + j*10e2; %impedancia de terminacion tubo abierto %Z = 100e-1 + j*100e-2; %Y = 1/Z; C = zeros(N,N); C(N,N) = 1; %C = ro*Y*A*C; %% incorporacion de la fuente en x = -L/2 %tanto para la matriz de masa como la mattriz de rigidez extraeremos las %matrices necesarias para este nuevo sistema M1 = M(2:N,2:N);%extrae de la matriz principal la submatriz delde la fila/columna 2 a la fila/columna N K1 = K(2:N,2:N);%extrae de la matriz principal la submatriz delde la fila/columna 2 a la fila/columna N C1 = C(2:N,2:N);%extrae de la matriz principal la submatriz delde la fila/columna 2 a la fila/columna N %% problema de valores y vectores propios % %sirve para determinar las frecuencias naturales y la distribucion de la % %presion sonora a dichas frecuencias naturales % %resolvemos la ecuacion % %[lam*M1 + K1]*phi = 0 % %M es matriz (N-1 X N-1) % %K es matriz (N-1 X N-1) % %phi es vector (N-1 X 1) % %lam es escalar lam = -omega^2 (1 X 1) % %fin de que la solucion sea no trivial, es decir distinta de cero % %obtendremos N soluciones para lam (es decir N frecuencias naturales) y N % %soluciones para phi(es decir N distribuciones de presion sonora) % % %eig soluciona el problema de valores propios % [PHI1, LAM1] = eig(K1,M1); % % %PHI es la matriz (N-1 X N-1)formada por los vectores propios % %PHI = %[phi1, phi2,..., phiN-1] los phi son vectores columna % %cada phi es (N-1 X 1) % % %LAM1 es la matriz diagonal formada por los valores propios % %LAM1 = diag(lam1, lam2, ..., lamN-1) % %LAM1 = diag(omega1^2, omega2^2, ..., omegaN-1^2) % % %frecuencias naturales aproximadas por FEM % %diag estrae la diagonal principal de una matriz % fFEM1 = real( sqrt (diag(LAM1) ) )/(2*pi); % % %normalizacion de vectores propios % mnorm = diag(transpose(PHI1)*M1*PHI1); % for n1 = 1:N-1 % PHI1(:,n1) = PHI1(:,n1)/sqrt( mnorm(n1) ); % end; %% resultados frecuencias % %resultados teoricos % fTEO1 = zeros(N-1,1); % for n1 = 1:N-1 % fTEO1(n1) = (2*n1-1)*c/(4*L); % end; % % %error en frecuencias % errof1 = zeros(N-1,1); % for n1 = 1:N-1 % errof1(n1) = 100*abs( (fTEO1(n1) - fFEM1(n1))/fTEO1(n1) ); % end; % % %grafico de frecuencias % figure(1) % plot(fTEO1,fFEM1,'-o'); % title('frecuencias teoricas vs frecuencias fem'); % xlabel('fTEO (Hz)'); % ylabel('fFEM (Hz)'); % grid on; % box on; % % %grafico de error % figure(2) % plot(errof1,'-o'); % title('error'); % xlabel('numero de modo'); % ylabel('error (%)'); % grid on; % box on; % % TablaError1 = [fTEO1 fFEM1 errof1]; %% resultados distribucion modal de presion sonora % %cordenadas % dx = L/Ne; %incremento o longitud de cada elemento % xcorde = (1:N-1)*dx - L/2; %nos da las coordenadas desde -L/2 a L/2 % % %grafico % figure(3) % % subplot(2,2,1) % plot( xcorde, PHI1(:,1)/max( abs( PHI1(:,1) ) ) ); % title('1er modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2 L/2 -1 1]); % grid on; % box on; % % subplot(2,2,2) % plot( xcorde, PHI1(:,2)/max( abs( PHI1(:,2) ) ) ); % title('2er modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2 L/2 -1 1]); % grid on; % box on; % % subplot(2,2,3) % plot( xcorde, PHI1(:,3)/max( abs( PHI1(:,3) ) ) ); % title('3er modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2 L/2 -1 1]); % grid on; % box on; % % subplot(2,2,4) % plot( xcorde, PHI1(:,4)/max( abs( PHI1(:,4) ) ) ); % title('4er modo FEM normalizado'); % xlabel('x (m)'); % ylabel('p (Pa)'); % axis([-L/2 L/2 -1 1]); % grid on; % box on; %% solución forzada con fuente en x = -L/2 %resolveremos el problema [-w^2*M1 + j*w*C1 + K1]P1 = F1 %esto consiste en resolver este sistema de ecuaciones para cada valor de %la frecuencia angular w = 2*pi*f, donde f = 1,2,...,1000 (Hz) es la frecuencia %la presión sonora en x = -L/2 es P(-L/2,t) = PL*exp(j*w*t) con PL = 1 (Pa) %armaremos los vectores de fueza con las primeras columnas de masa y %rigidez originales, es deccir antes de colocar las condiciones de contorno PL = 1; Pm = M(2:N,1); Pk = K(2:N,1); Pc = C(2:N,1); %frecuencias fmin = 1; fmax = 1000; finc = 1; f = fmin:finc:fmax; %frecuencia w = 2*pi*f; %frecuencia angular NF = length(f); %longitud vector frecuencia %prelocalizacion de los resultados de la presion sonora %se guardan los resultados en una matriz de N-1 filas que indican la %posicion en el eje "x" y NF columnas que son para las frecuencias P1 = zeros(N-1,NF); %la prelocalización nos ayuda con la velocidad de ejecución %ciclo principal de calculo for nf = 1:NF %numero de onda k = w(nf)/c; %impedancia de radiacion ka2 = 2*k*rd; RL = 1 - 2*besselj(1,ka2)/ka2; XL = 2*struve(1,ka2)/ka2; ZL = ro*c*( RL + j*XL); %matriz de amortiguamiento C1(N-1,N-1) = ro*A/ZL; %vector de fuerza F1 = w(nf)^2*PL*Pm - j*w(nf)*PL*Pc - PL*Pk; %vector de fuerza dependiente de la frecuencia %resolucion P1(:,nf) = (-w(nf)^2*M1 + j*w(nf)*C1+ K1)\F1; %A\b es lo mismo que inv(A)*b pero mas exacta end; %despues debo incorporar para todas las frecuencias la condición de %contorno p(-L/2) = PL Pini(1,1:NF) = PL; %fila con valores PL para todas las frecuencias P = [Pini;P1]; %concatemanos la fila Pini a la matriz P1 %% resultados %compararemos resultados de la distribucion de la presion sonora %para una frecuencia especifica que no sea de resonancia f = 1000 (Hz) %en la longitud del tubo %posicion x a partir de x = 0 dx = L/Ne; %incremento o longitud de cada elemento xcorde = (1:N-1)*dx - L/2; %nos da las coordenadas desde -L/2 a L/2 x = xcorde + L/2; x = [0,x]; NX = length(x); %numero de onda k = 2*pi*f(500)/c; %resultado teórico tubo cerrado para impedancia alta for nx = 1:NX pteo(nx) = PL*cos(k*(L - x(nx)))/cos(k*L); end; %resultado teórico tubo abierto para impedancia baja % for nx = 1:NX % pteo(nx) = PL*sin(k*(L - x(nx)))/sin(k*L); % end; figure(4) plot(x,pteo,'b',x,P(:,500),'r -o') title('comparacion resultados teoricos vs FEM'); xlabel('x (m)'); ylabel('p (Pa)'); legend('pteo','pFEM'); grid on; box on; %error cuadratico erro2prom = sum( abs((pteo - transpose(P(:,500)) )).^2 ); %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(P),'FaceColor','interp',... 'EdgeColor','none',... 'FaceLighting','phong') title('Modulo Presion Sonora - Tubo Cerrado Impedancia - FEM'); xlabel('f (Hz)'); ylabel('x (m)'); zlabel('|P| (Pa)') axis([0,f(NF),0,L]); grid on; box on; colormap jet colorbar; figure(6) surf(Fmesh,Xmesh,20*log10(abs(P)/sqrt(2)/2e-5),'FaceColor','interp',... 'EdgeColor','none',... 'FaceLighting','phong') title('Nivel de Presion Sonora - Tubo Cerrado Impedancia - FEM'); xlabel('f (Hz)'); ylabel('x (m)'); zlabel('|P| (Pa)') axis([0,f(NF),0,L]); grid on; box on; colormap jet colorbar; figure(7) plot(f,20*log10(abs(P(N,:))/sqrt(2)/2e-5)); title('Solucion FEM '); title('respuesta de frecuencia al final del tubo'); xlabel('f (Hz)'); ylabel('Lp (dB)'); legend('Lp'); grid on; box on;