%% prog fem012 tubo fuente en un extremo impedancoia radiacion %este programa no calcula los valores propios de un tubo de %longitud L con fuente en un extremo con impedancia de radiación %area de la seccion transversal constante A %el numero de elementos es variable %se utilizara matrices elementales clear all; close all; clc; %% inicializacion 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 = 1e-4; %area seccion transversal tubo cuadrado 0.01 X 0.01 m2 c = 340.0; %velocidad del sonido ro = 1.12; %densidad aire B = c^2*ro; %rigidez de compresibilidad adiabática de volumen rd = sqrt(A/pi); %radio de tubo; rd1 = sqrt(A/pi)/5; %radio de otro hoyo en la mitad del tubo A1 = pi*rd1^2; %area seccion transversal apertura en la mitad %% 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 elemental y global % la matriz de masa es normalizada a la velocidad del sonido %matriz de masa elemental Me = [2 1;1 2]; Me = (A*ax)/(3*c^2)*Me; %inicializacion M = zeros(Nt,Nt); %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 ); end; end; end; end; end; end; %% matriz de rigidez %matriz de rigidez elemental Ke = [1 -1;-1 1]; Ke = A/(2*ax)*Ke; %inicializacion K = zeros(Nt,Nt); %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; K(NN,MM) = K(NN,MM) + Ke( fila, colu ); end; end; end; end; end; end; %% 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(Nt,Nt) = 1; %% incorporacion de la fuente en x = -L/2 %tanto para la matriz de masa como la matriz de rigidez extraeremos las %matrices necesarias para este nuevo sistema M1 = M(2:Nt,2:Nt);%extrae de la matriz principal la submatriz desde la fila/columna 2 a la fila/columna N K1 = K(2:Nt,2:Nt);%extrae de la matriz principal la submatriz desde la fila/columna 2 a la fila/columna N C1 = C(2:Nt,2:Nt);%extrae de la matriz principal la submatriz desde la fila/columna 2 a la fila/columna N %% problema de valores propios % % [PHI1,LAM1] = eig(K1,M1); % % %frecuencias aproximadas % fFEM1 = real( sqrt( diag(LAM1) ) )/(2*pi); %formato compactor % % %normalizacion de los vectores propios % mnorm = diag( transpose(PHI1)*M1*PHI1 ); % % for n1 = 1:Nt-1 % PHI1(:,n1) = PHI1(:,n1)/mnorm(n1); % end; % % %frecuencias teoricas % for n2 = 1:Nt-1 % fTEO(n2,1) = (2*n2-1)*c/(4*L); % end; % % %para visualizar el error % Nerr = Nt-1; % % %frecuencias teoricas y FEM % figure(1) % plot(fTEO(1:Nerr),fFEM1(1:Nerr),'-o'); % title('Frecuencias Teoricas vs Frecuencias FEM'); % xlabel('f_{TEO} (Hz)'); % ylabel('f_{FEM} (Hz)'); % grid on; % box on; % grid on; % % %% calculo del error % % %error % errfrec(1:Nt-1,1) = 0; % for n1 = 1:Nt-1 % errfrec(n1) = 100*abs(fTEO(n1) - fFEM1(n1) )/fTEO(n1); % end; % % % %tabla % TablEerror = [fTEO fFEM1 errfrec]; % % %eror % figure(2) % plot(errfrec(1:Nerr),'-o'); % title('Error Frecuencias Teoricas vs Frecuencias FEM'); % xlabel('N de Modo'); % ylabel('error [%]'); % grid on; % box on; % grid on; % % % %% modos normales % % pini = zeros(1,Nt-1); % PHI1 = [pini;PHI1]; % % %modos normales % figure(3) % subplot(2,2,1) % plot( xcorde, PHI1(:,2)/max(PHI1(:,2)) ) % title('1er Modo FEM'); % grid on; % box on; % grid on; % axis([-L/2 L/2 -1 1]); % % subplot(2,2,2) % plot( xcorde, PHI1(:,3)/max(PHI1(:,3)) ) % title('2do Modo FEM'); % grid on; % box on; % grid on; % axis([-L/2 L/2 -1 1]); % % subplot(2,2,3) % plot( xcorde, PHI1(:,4)/max(PHI1(:,4)) ) % title('3er Modo FEM'); % grid on; % box on; % grid on; % axis([-L/2 L/2 -1 1]); % % subplot(2,2,4) % plot( xcorde, PHI1(:,5)/max(PHI1(:,5)) ) % title('4to Modo FEM'); % grid on; % box on; % grid on; % axis([-L/2 L/2 -1 1]); %% solución forzada con fuente en x = -L/2 %% 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 = 1e-2; 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 apertura final ka2 = 2*k*rd; RL = 1 - 2*besselj(1,ka2)/ka2; XL = 2*struve(1,ka2)/ka2; ZL = ro*c*( RL + j*XL); %impedancia de radiacion apertura mitad ka21 = 2*k*rd1; RL1 = 1 - 2*besselj(1,ka21)/ka21; XL1 = 2*struve(1,ka21)/ka21; ZL1 = ro*c*( RL1 + j*XL1); %matriz de amortiguamiento C1(N-1,N-1) = ro*A/ZL; C1( round( (N-1)/2 ) , round( (N-1)/2 ) ) = ro*A1/ZL1; %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;