4.6 Computer Program |
105 |
Complex specific acoustic impedance for plane waves associated with e2jωt
(e jωt) |
p (x, t); u (x, t) |
|
|
|
|
|
|
|
|
|
|
|
|
z (x) ¼ p /u |
|
|
|
||||||||
p +p |
p : [(A+c |
|
jA+s)e jkx |
+ (A |
|
c |
|
|
jA s)e jkx]e jωt |
|
A c |
|
jA s ejkx |
A c |
|
jA s e jkx |
|||||||||
þ |
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
ρoc |
ð þ |
þ Þ þð |
Þ |
|||||||
|
|
1 |
|
|
|
jkx |
|
|
jkx |
jωt |
A c |
|
jA s ejkx |
A c |
|
jA s e jkx |
|||||||||
u |
+u |
u : |
[(A+c |
|
jA+s)e |
|
|
(A |
|
c |
|
jA |
|
s)e |
]e |
|
ð þ |
|
þ Þ ð |
|
Þ |
||||
|
|
|
|||||||||||||||||||||||
þ |
2 |
|
ρo c |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
Terms associated with e jωt are the conjugates of the terms associated with e jωt. Usually, the analysis of sound transmission in pipes only needs to compute for one part associated with e jωt. The results of p* and u* associated with e jωt are just
the conjugates of the resulting p and u using the complex z associated with e jωt. The final real solutions p and u will be the summation of the pair of the complex
solution.
4.6Computer Program
Plotting Traveling and Standing Waves
A standing wave, ps(t, x), is the addition of a forward wave, p+(t, x), and a backward wave, p (t, x), as shown below:
p+(t, x) ¼ Ac cos (ωt kx) As sin (ωt kx) ¼ A cos (ωt kx + ϕ)
p (t, x) ¼ Ac cos (ωt + kx) As sin (ωt + kx) ¼ A cos (ωt + kx + ϕ)
ps(t, x) ¼ p+ + p
¼ 2A cos (ωt + ϕ) cos (kx)
(Form 1:RIP) (Form 2:REP) (Form 1:RIP) (Form 2:REP)
See Exercise 3.4
Use MATLAB to plot the pressure and velocity of a forward traveling wave (Part A), a backward traveling wave (Part B), and a standing wave (Part C) of three wavelengths using the following parameters:
Ac ¼ 0 ½m&; As ¼ 2 ½m&; ω ¼ 2π |
s |
; k ¼ 2π |
m |
||
|
|
rad |
|
|
rad |
Plot eight subplots to capture the wave motion with a time increment of T divided by eight. For calculating the velocity, assume ρoc ¼ 1.s
Part A:
•Study and understand the provided MATLAB code.
•Complete the function “getForwardWave ” which calculates the pressure and velocity from A, w, k, Time, and XDir.
•Run the code to plot the following eight subplots representing the forward wave motion during a period T:
pþðt, xÞ ¼ A cos ðωt kx þ ϕÞ
106 |
4 Acoustic Intensity and Specific Acoustic Impedance |
Sound Pressure |
Molecular Flow Velocity |
Part B:
•Complete the function “getBackwardWave” which calculates the pressure and velocity from A, w, k, Time, and XDir.
•Run the code to plot the following eight subplots representing the backward wave motion during a period T:
p ðt, xÞ ¼ A cos ðωt þ kx þ ϕÞ
Sound Pressure Molecular Flow Velocity
4.6 Computer Program |
107 |
Part C:
•Complete the function “getStandingWave ” which calculates the pressure and velocity from A, w, k, Time, and XDir.
•Run the code to plot the following eight subplots representing the standing wave motion during a period T:
psðt, xÞ ¼ pþ þ p
|
¼ 2A cos ðωt þ ϕÞ cos ðkxÞ |
Sound Pressure |
Molecular Flow Velocity |
MATLAB Code:
function ANC_PRJ021_TravelingAndStandingWaves |
|
|
clear all % removes all variables |
|
|
close all % deletes all figures |
|
|
%---------------------------------------- |
–-------- |
–% |
% Section 1: Define the Variables and Parameters |
|
|
%------------------------------------------------- |
|
–% |
%A standing wave,ps(t,x),is the addition of a forward wave, p+(t,x),
%and a backward wave, p-(t,x) as shown below
% p+(t,x)=Ac cos(wt-kx)+As sin(wt-kx) |
(Form 1) |
|
% |
=A cos(wt-kx+phi) |
(Form 2) |
% p-(t,x)=Ac cos(wt+kx)+As sin(wt+kx) |
(Form 1) |
|
% |
=A cos(wt+kx+phi) |
(Form 2) |
% ps(t,x)=P+(t,x)+p-(t,x) |
|
|
108 |
4 Acoustic Intensity and Specific Acoustic Impedance |
|
% |
=A cos(wt-kx+phi)+A cos(wt+kx+phi) |
(Form 2) |
% |
=2A cos(wt+phi)*cos(kx) --- see the Homework Eroblem 3.4 |
|
%note that the p+ and p- have the same amplitude A
%define the parameters of the wave function iCase=1;
if iCase == 1 % phase phi = 0
Ac=0;% the coefficient of the cosine function in Form 1 As=2;% the coefficient of the sine funciton in Form 1 elseif iCase == 2
Ac=3^(1/2); % the coefficient of the cosine function in Form 1 As=-1; % the coefficient of the sine funciton in Form 1
end
w=2*pi/1; |
% the angular frequency [rad/sec] in Form 1 and Form 2 |
|
k=2*pi/1; |
% the wave number [rad/m] in Form 1 and Form 2 |
|
% define the number of times to take a snapshot |
||
nTimePerPeriod=8; |
% the number of data per period (T) |
|
%plot only one period since it repeats every period
%define the parameters for plotting [space related]
nDataPerWaveLength=8; |
% the number of data per period (T) |
||
nWaveLength=3; |
% the number of periods for plotting |
||
%----------------------------------------------------- |
|
|
–% |
% Section 2: Setup |
|
|
|
%---------------------------------------------------------- |
|
|
–% |
% get the period (T) and a time array |
|
||
T =2*pi/w; |
|
% the period of the wave [sec] |
|
TimeIncrement=T/nTimePerPeriod; |
% [sec] |
||
nTime=nTimePerPeriod; |
% the number of time data |
||
TimeArray=(0:1:nTime-1)'*TimeIncrement; % [sec]
%get the data array of the location in the X-direction L=2*pi/k; % [m] , note that k=2*pi/L XDataIncrement=L/nDataPerWaveLength; % [m]
nXData=nDataPerWaveLength*nWaveLength; % the number of space data XDataArray=(0:1:nXData-1)'*XDataIncrement; % [m]
%get the data array of the location in the Y-direction
%(Y-direction is a dummy direction and does not exist) YDataIncrementDummy=XDataIncrement*0.15; % [m] a dummy value nYDataDummy=5; % a dummy number(to adjust the arrow size) YDataArrayDummy=(0:1:nYDataDummy-1)'*YDataIncrementDummy; % [m]
%get A and phi in Form 2 from From Ac and As in Form 1
[A,phi]=get_A_and_phi(Ac, As); % the same as PROJECT 1.9.1 |
|
%-------------------------------------------------------------- |
–% |
% Section 3: Calculation |
|
%-------------------------------------------------------------- |
–% |
% Index of Cases |
|
% iCaseWaveDirection = 1; % Forward Traveling Wave |
|
% iCaseWaveDirection = 2; % Backward Traveling Wave |
|
% iCaseWaveDirection = 3; % Standing Wave |
|
for iCaseWaveDirection=1:1:3 |
|
% calculate the cosine function |
|
for kTime=1:1:nTimePerPeriod |
|
for kXData=1:1:nXData |
|
% extract the scalars from the arrays |
|
% scalars are used temporally in the for-loop |
|
4.6 Computer Program |
109 |
TimeTemp=TimeArray(kTime,1);
XDataTemp=XDataArray(kXData,1);
if iCaseWaveDirection == 1 % forward wave
% get the pressure from A and phi |
(Form 2) |
[pressTemp,velocTemp]=... |
|
getForwardWave(A,w,k,TimeTemp,XDataTemp,phi); |
|
% put the pressure and velocity in the matrices pressColMatrix(1,kXData)=pressTemp; % pressure veloXColMatrix(1,kXData)=velocTemp; % velocity veloYColMatrix(1,kXData)=0; % dummy value % adjust the arrow size in the figures
arrow_scale_factor=1/20; % adjust the arrows size elseif iCaseWaveDirection == 2 % backward wave
[pressTemp,velocTemp]=...
getBackwardWave(A,w,k,TimeTemp,XDataTemp,phi); % put the pressure and velocity in the matrices pressColMatrix(1,kXData)=pressTemp; % pressure veloXColMatrix(1,kXData)=velocTemp; % velocity
veloYColMatrix(1,kXData)=0; |
% dummy value |
% adjust the arrow size in figure |
|
arrow_scale_factor=1/20; |
|
elseif iCaseWaveDirection == 3 % standing wave |
|
[pressTemp,velocTemp]=... |
|
getStandingWave(A,w,k,TimeTemp,XDataTemp,phi); |
|
pressColMatrix(1,kXData)=pressTemp; % pressure |
|
veloXColMatrix(1,kXData)=velocTemp; % velocity |
|
veloYColMatrix(1,kXData)=0; |
% dummy value |
% adjust the arrow size in the figures |
|
arrow_scale_factor=1/40; |
|
end |
|
end |
|
%------------------------------------------------------------ |
% |
% Section 4: Plotting |
|
%------------------------------------------------------------ |
% |
%the meshed data of the pressure for the surface plots onesRowMatrix=ones(nYDataDummy,1); pressGrid=onesRowMatrix*pressColMatrix; veloXGrid=onesRowMatrix*veloXColMatrix; veloYGrid=onesRowMatrix*veloYColMatrix;
%figures setup
figure(210+iCaseWaveDirection); set(gcf,'Position',[100 100 1000 600]) subplot(8,2,kTime*2-1) plot(XDataArray,pressColMatrix,'lineWidth',2) xlim([0 3])
xticks([0 0.5 1 1.5 2 2.5 3]) xticklabels({'0','\lambda/2','\lambda','3\lambda/2',...
'2\lambda','5\lambda/2','3\lambda'}) ylim([-4 4])
yticks([-4 0 4]) str=strcat('(',num2str(kTime-1),'/8)T'); ylabel(str)
% plot the pressure colormap
110 |
4 Acoustic Intensity and Specific Acoustic Impedance |
|
subplot(8,2,kTime*2-0) |
|
pcolor(XDataArray,YDataArrayDummy,pressGrid); shading interp; |
|
xlim([0 3]) |
|
xticks([0 0.5 1 1.5 2 2.5 3]) |
|
xticklabels({'0','\lambda/2','\lambda','3\lambda/2',... |
|
'2\lambda','5\lambda/2','3\lambda'}) |
|
set(gca,'YTick',[]) |
|
cb=colorbar('location','eastoutside'); |
|
set(cb,'position',[.92 .125 .01 .795]) |
|
caxis([-4 4]) |
|
hold on |
|
% plot the velocity vectors |
|
subplot(8,2,kTime*2-0) |
|
h=quiver(XDataArray,YDataArrayDummy,... |
|
veloXGrid.*arrow_scale_factor,... |
|
veloYGrid.*arrow_scale_factor); |
|
xlim([0 3]) |
|
ylim([0 max(max(YDataArrayDummy))]) |
|
set(gca,'YTick',[]) |
|
set(h,'AutoScale','off') |
|
set(h,'LineWidth',1) |
|
set(h,'color','w') |
|
hold off |
end |
|
% save the figure |
|
if iCaseWaveDirection == 1 % forward wave |
|
saveas(gcf,strcat('Figure_ANC_PRJ021_ForwardWaves'),'emf') |
|
elseif iCaseWaveDirection == 2 % backward wave |
|
saveas(gcf,strcat('Figure_ANC_PRJ021_BackwardWaves'),'emf') |
|
elseif iCaseWaveDirection == 3 % standing wave |
|
saveas(gcf,strcat('Figure_ANC_PRJ021_StandingWaves'),'emf') |
|
end |
|
end |
|
end |
|
%-------------------------------------------------------------- |
–% |
% Section 5: Functions |
|
%-------------------------------------------------------------- |
–% |
function [pressure,velocity]=getForwardWave(A,w,k,Time,XDir,phi)
% p+(t,x)=Ac cos(wt-kx)+As sin(wt-kx) |
(Form 1) |
|
% |
=A cos(wt-kx+phi) |
(Form 2) |
%modify the following two lines for the pressure and velocity
%assume loc=0
pressure=A*cos(w*Time-k*XDir+phi); % pressure of a forward wave
velocity=pressure; |
% velocity of a forward wave |
end |
|
function [pressure,velocity]=getBackwardWave(A,w,k,Time,XDir,phi)
% p-(t,x)=Ac cos(wt+kx)+As sin(wt+kx) |
(Form 1) |
|
% |
=A cos(wt+kx+phi) |
(Form 2) |
%modify the following two lines for the pressure and velocity
%assume loc=0
pressure=0; % pressure of a backward wave velocity=0; % velocity of a backward wave end