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 

   c_right = 2*c; % speed to the right of the disc
   for i=1:N
      if X(i) > D
         C(i) = c_right;
   % 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);

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

1 Answer 1


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.

  • $\begingroup$ thank you. I am afraid I am still a little confused. Given that the density is not independent of x at the boundary does this mean that the wave equation does not apply at the boundary? Also what does the variable "v (or u)" refer to? Would it be the velocity of the particle and if so why should this be constant across the interface. Thanks again. I really appreciate your help $\endgroup$
    – Jacob_K
    Commented Apr 12, 2020 at 21:06
  • $\begingroup$ 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! $\endgroup$ Commented Apr 12, 2020 at 22:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.