%% Prog001FEMTuboCerradoCerradoSinMontaje clearvars; close all; clc; %este es un programa de elementos finitos que representa un tubo cerrado %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 este caso comenzaremos con 4 elementos 5 nodos %% datos inciales L = 1.0; %longitud del tubo Ne = 4; %numero de elementos N = Ne + 1; %numero de nodos grados de libertad DOF ax = L/(2*Ne); %semilongitud de cada 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 %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 %% matriz de masa 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]; M = (A*ax)/(3*c^2)*M; %% matriz de rigidez 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]; K = A/(2*ax)*K; %% valores y vectores propios %la función eig() determina los valores y vectores propios de la matriz [PHI,LAM] = eig(K,M); %frecuencias naturales elementos finitos FEM f_FEM = diag( real(sqrt(LAM)))/(2*pi); %frecuencias naturales teóricas TEO f_TEO = zeros(N,1); for n = 1:N f_TEO(n) = (n-1)*c/(2*L); end; %error err = 100*abs(f_FEM - f_TEO)./f_TEO; %tabla f_FEM_f_TEO_err = [f_FEM, f_TEO, err]; %% figuras frecuencias valores propios figure(1) plot(f_TEO,f_FEM,'-o'); title('frecuencias f_{TEO} vs f_{FEM}'); xlabel('f_TEO (Hz)'); ylabel('f_TEO (Hz)'); grid on; box on; figure(2) plot(err,'-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, PHI(:,1)/max( abs(PHI(:,1) ) ) ) ; title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(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, PHI(:,2)/max( abs(PHI(:,2) ) ) ) ; title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(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, PHI(:,3)/max( abs(PHI(:,3) ) ) ) ; title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(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, PHI(:,4)/max( abs(PHI(:,4) ) ) ) ; title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(4))] ); xlabel('x (m)'); ylabel('p (Pa)'); axis([-L/2,L/2,-1.5,1.5]); grid on; box on;