# Simulating pressure waves at an impedance boundary

I am trying to simulate pressure waves crossing a boundary from one medium to another (e.g., water to air) in Matlab. The code that I have got so far, which is largely taken from Wikipedia on Partial Transmission only changes the speed as the wave crosses the boundary. The code works using a recursive form or the wave equation found using the method of central differences.

I am confused about how to incorporate the effect of density on the impedance of the medium. From my understanding:

$$Z=\rho C$$ where $$Z$$ is the impedance, $$\rho$$ is the density, $$C$$ is the speed of sound

At the boundary, the proportion of the waves' amplitude that is reflected and transmitted is given by the reflection and transmission coefficients which are: $$R=\frac{Z_1-Z_2}{Z_1+Z_2},\\ T=R=\frac{2Z_1}{Z_1+Z_2}$$

At the moment the simulation will correctly calculate the reflected and transmitted waves amplitudes for two mediums with the same density but different wave speeds; however, I have no way of varying the density of the medium.

Seeing as neither density or impedance are relevant in the wave equation I am confused about how they end up determining the reflection and transmission coefficients and how to implements them into my simulation.

Below is the relevant section of my code:

   % length of the string and the grid
L = 5;
N = 151;
X=linspace(0, L, N);

h = X(2)-X(1); % space grid size
c = 0.01; % speed of the wave for visualisation
tau = 0.25*h/c; % time grid size

% form a medium with a discontinuous wave speed
C = 0*X+c;  %this has formed a vector the same dimension as X with every entry =c

D=0.5*L;
c_right = 2*c; % speed to the right of the disc
for i=1:N
if X(i) > D
C(i) = c_right;
end
end
% Now C = c fo x < D, and C=c_right for x > D

K = 10; % steepness of the bump
S = 0; % shift the wave
f=inline('exp(-K*(x-S).^2)', 'x', 'S', 'K'); % a gaussian as an initial wave
df=inline('-2*K*(x-S).*exp(-K*(x-S).^2)', 'x', 'S', 'K'); % derivative of f

% wave at time 0 and tau
U0 = 0*f(X, S, K);
U1 = U0 - 2*tau*c*df(X, S, K);

U = 0*U0; % current U

% plot between Start and End
Start=0; End=1500;

for j=1:End

%  fixed end points
U(1)=0; U(N)=0;

% finite difference discretization in time
for i=2:(N-1)
%this is the wave equation written algebraically with second
%order central difference theorem and then rearranged for U(i)
U(i) = (C(i)*tau/h)^2*(U1(i+1)-2*U1(i)+U1(i-1)) + 2*U1(i) - U0(i);
end

% update info, for the next iteration
U0 = U1; U1 = U;


Essentially going one step back in the derivation of the wave equation you take the divergence of the (simplified) Euler equation $$\mathbf{\nabla}.(\rho_{0} \frac{d\mathbf{v}}{dt}=-\nabla p$$) and next, you let $$\rho_{0}$$ out of the operator, assuming it is a constant. Now at the inferface between two media with $$\rho_{01} \neq \rho_{02}$$, $$\nabla \rho_{0}$$ is not null anymore, so this simplification is not true. That's the physics part. So you need to implement a specific boundary condition at the interface between two such media in your model to ensure both $$p$$ and the normal component of particule velocity $$v$$ (or displacement $$u$$) are continuous accross the interface.
• Sorry I corrected my notations, $v$ now consistently represents particule velocity, $u$ is the displacement. Also the right boundary condition is that the normal component of velocity( or displacement) should be continuous across the interface (look at how the reflection and transmission coefficients are derived). The normal component may be obtained form Euler equation ($\rho dv_{x}/dt= -dp/dx$) if boundary normal is along $x$, for instance). This might not be a simple condition to implement in your numerical scheme, but I cannot be of any further help, I'm afraid! – user8736288 Apr 12 '20 at 22:08