Engineering Acoustics/Car Mufflers:Animation:code

% 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