411 REFERENCE · Wikibooks
Engineering Acoustics/Car Mufflers:Animation:code
Reference material by Wikibooks contributors. This retains its source identity and is not 411 news reporting.
% Duct with enlargement - last modified 11/03/99
% Reproduit avec l'aimable permission de Prof. L. Mongeau
clear all;
clc;
close all;
Nsteps=30;
f=86*2; % frequency
rho=1.15;
c=344.;
rhoc=rho*c;
l=1; %length of duct
a1=0.0762/2;
S1=pi*a1^2;
a2=0.1524/2;
S2=pi*a2^2;
m=S2/S1;
w=2*pi*f;
k=w/c;
t=linspace(2*pi/Nsteps/w,2*pi/w,Nsteps); %Nsteps time steps per period
alpha1=2.93e-5/a1.*sqrt(f);
alpha2=2.93e-5/a2.*sqrt(f);
khat=k-i*alpha1;
kl=khat*l;
U0=0.01;
R=(1-m)/(m+1);
A=rho*c*U0./(exp(i*kl)-R*exp(-i*kl));
B=R*A;
C=A+B;
x=linspace(0,2,100); % computational domain
y=linspace(-0.2,0.2,32);
for m=1:50
for k=8:24
pt(k,m)=A*exp(i*khat*(l-x(m)))-B*exp(-i*khat*(l-x(m)));
end
end
for m=51:100
for k=1:32
pt(k,m)=C*exp(i*khat*(l-x(m)));
end
end
pmag=abs(pt);
pphase=angle(pt);
clim=max(max(pmag));
V=[-clim clim]; % sets color scheme
pt=pmag.*exp(i*pphase); % red: maximum positive
% blue: minimum negative pressure
figure(1)
pt2=real(pt);
H=pcolor(x,y,pt2);
axis equal % to eliminate distortion
axis([0 2 -.2 .2])
M=moviein(Nsteps); % for Matlab movie animation
for k=1:Nsteps
pt2=real(pt*exp(i*w*t(k)));
H=pcolor(x,y,pt2);
axis equal
xlabel('x');
ylabel('y');
title('pressure')
shading interp;
caxis(V);
axis([0 2 -.2 .2])
M(:,k)=getframe;
% to store the images for web posting (bitmaps)
% filnme=strcat('spher',int2str(k));
% eval(['print -dbitmap ' filnme ' -f']);
end
movie(M,20)
Back to Car Mufflers