NNESMO Practical 3: The wave equation

Prof Stefan Nielsen

Department of Earth Sciences

Durham University

Introduction
Part 1: Harmonic oscillator
Part 2: 1D wave on a spring
Part3: 2D wave
Part 4 (advanced): higher order interpolation
Further notes on higher order interpolation:

Introduction

In the three practicals below, increasing complexity is incorporated by adding an extra dimension each time. The problem solved is very similar in the three cases, involving an oscillation in time, which propagates in the 1D and 2D cases.

We start (part 1) from zero dimensions — harmonic stationary oscillation, i.e. the behaviour of a simple spring and mass system; no propagation is involved in this case. This problem has well known and simple mathematical solutions, so the numerical implementation only serves as a preliminary exercise (warming up) and as an illustration of the concept of stability, discretisation and sampling rate.

We then will move to a propagating oscillation (part 2) along a single direction (i.e. 1D propagation similar to waves along a string). In this case we will start to build a simple grid to locate the sample points in both space and time. The concept of stability is further explored, with the introduction of the CFL criterion (see lecture slides).

Finally, we will implement the wave propagation in 2D space plus time (propagating wave across a section of the Earth, a problem similar to that of a vibrating membrane).

In all cases, a single displacement is considered, and the problem is exactly equivalent to that of a pressure wave (acoustic). Note that in the 2D case, this “acoustic” problem can also represent the SH mode of propagation of elastic waves in the Earth (Fig 1). In SH mode (Shear waves, Horizontal motion) the particles move in the horizontal direction (call it ) perpendicular to the vertical section of the model (call it ). In this case there are two shear stresses (equivalent to the two gradients of pressure in the and in the directions, for the case of the acoustic wave problem). For simplicity, we may call or the stress arising from the horizontal (along ) motion derivative, and or that arising from the vertical (along ) motion derivative.

Part 1: Harmonic oscillator

3.1) Write the equation for a harmonic oscillator in the case of a mass and spring (i.e., acceleration times mass equal displacement times stiffness). Then write a basic analytical solution using a periodic function with angular frequency . Replace it in the equation and find the period of oscillation in terms of mass and stiffness.

3.2) Now re-write the oscillator equation as a system of two coupled equations, one for velocity time derivative, the other for the force time derivative (see slides).

3.3) Replace the differentials in the two equation system above with their discrete version (time update), for a time stepping . Build a time grid with the alternating velocity and force values.

3.4) Program the equations into matlab:

First you will have to initialise the parameters (mass, stiffness, time step…) then the variables (force, velocity…). For the variables, choose adequate initial conditions — if the mass is initially at rest and at equilibrium, it will not move! — so you need to choose an initial condition where either the mass has non-zero velocity, or the mass is off the equilibrium position ( ) so that it will accelerate.

The next thing to go in the matlab code is the time iteration loop, with the equations for force and velocity updating. The duration should be at least a few periods of the oscillation (10-20) so you can see what happens. You can update velocity first, then the force, or vice-versa, it does not really matter.

Plot the result. Also compute the basic periodic analytical solution at the end of the code and plot it against the numerical solution.

Try different time steps. What happens if is of the same order as the oscillation period or larger? What is the limiting value for a stable solution? What is the limiting value for a smooth solution?

ADVANCED: try to add viscous damping, that is, a force that opposes motion and is proportional to velocity. At which points of the grid (i.e. force points or velocity points) should that force be computed in principle? At which points should it be applied? How do you integrate it properly?

Throughout these practicals we can use the same notation and the same parameters for simplicity. For further development into waves (with application to seismology), I suggest we start to use mass and stiffness in this problem which may correspond to the mass and the elastic stiffness of a cubic meter of rock mass, i.e., GPa (i.e., a 1m spring of stiffness N) and kg/m3 (i.e., a cube of 1m side weighing 2.5 tons).

Template

Write the code from scratch or, if needed, use the following template and fill the ?? gaps.

% Initialization of parameters:
clear;
mu=1e10; rho=2500.; % stiffness and mass values
eta=1E6; % viscosity (optional for damping)
omega= ?? ; % angular frequency of the oscillator is sqrt(stiffness/mass)
T= ??; % period of the oscillator, by definition of period vs. frequency.
dt=T/(3*pi); % time stepping, a fraction of the period
nnT=10; %this defines how many periods the computation should last
nt=round(nnT*T/dt); % this computes the number of time iterations nt
vdf=dt*eta/(2*rho); % this is to pre-compute the viscous damping factor
%
%
%
% Initialization of the variables:
v= ?? ; s= ?? ; % initialise scalar variables (velocity & stress)
u(1:nt)= ?? ; % initialise vector variables (displacement)
%
%
% Start time loop
for t=1:nt
 % write the velocity update
 ??
 % write the stress update
 ??
 % update displacement and save for plotting
 ??
end
% plot the numerical solution:
figure (1);
??
% compute then plot the analytical solution:
??
% add graphic cosmetics:
title(['steps per period:' num2str(2*pi/(omega*dt))]);
set(gca,'FontSize',12,...
 'TickDir','out',...
 'TickLength',[.02,.02] )%, 'Ylim',[-1E-3 , 1E-3]);

Part 2: 1D wave on a spring

3.5) Write the 1D equation introducing rigidity and density , as shown in the lecture slides for the “spring” model. Working with the particular solution try to find the velocity of the waves in terms of by using , and the dispersion relation .

3.6) Now re-write the above wave equation as a system of two coupled equations, one for velocity time derivative, the other for the force time derivative (see slides).

3.7 Define boundary conditions – a 1D medium has two (one a the beginning, and one at the end, ). You may define one as a “fixed” or rigid boundary, and the other as a “free” or zero stress boundary.

3.8) Replace the differentials in the two equation system above by their discrete version (time update), for a time stepping . Build a time and space grid where the velocity values and the stress values are defined on different nodes so that derivatives are centred; this should create a particular form of staggered grid with “empty” or dummy nodes, and other nodes with values defined. See lecture slides if you experience difficulty, anxiety, unpleasant bowel movements or uncontrollable laughter. If all else fails, ask help from the assistants or lecturer. Do not, under any circumstance, throw the computer screen out of the window.

3.9) If the spatial stepping is m, what typical wavelength can you propagate with such sampling? What is the corresponding frequency (use the dispersion relation)?

3.10) Program into matlab. Use the same sequence of steps as in the oscillator, namely:

Start by defining all parameters and initialising the variables. You may use the following dimensions / parameters:

– Stiffness and density values : ; kg/m3.

– m

– Apply the CFL criterion to determine the maximum value of time stepping as a function of and (use the formula for as a function of density and rigidity as in point 1.1 above).

– Space and time number of points: ; .

Proceed to the time loop and apply the update scheme to both velocity and stress. You may use the implicit shifting (a) or the diff (b) function to calulate derivatives of arrays:

(a) u(2:nx)-u(1:nx-1) which is equivalent to u(i+1)-u(i) at all i.

(b) diff(u) which is equivalent to the above (see pre-course preparation notes).

Within the time loop you will also have to integrate the boundary conditions; use rigid B.C. ( ) at and free B.C. ( ) at .

Design a source time-function whose dominant frequency is compatible with the spatial grid. You may use a Gaussian in time, with a characteristic duration , e.g., . Assume that is the dominant period of the source, . Using the dispersion relation , which can be written also as , set so that the wavelength is at least ten times . The source perturbation may be added to either the velocity or the stress component in one point of the grid. Take care that the gaussian source is offset by a sufficient time interval so that its initial portion is not truncated.

Plot the results; you may dynamically plot the motion velocity, during the time loop, by updating the figure every 50 time iterations (use a test on a incremented integer variable) to visualise the propagating wave(s). Also plot the wave in an (x, t) graphic after the end of the time loop.

What happens to the wave at the rigid boundary, and at the free boundary?

Now change the properties of one part of the string only. To do that you will need to define and as a function of position in the model, so in the code they will be a function of the grid index and they have to be initialised as vectors instead of scalars. The model should reflect the following sequence of properties:

– Stiffness and density values for : ; kg/m3.

– Stiffness and density values for the : ; kg/m3.

– What happens at the interface between the two different string portions?

Template

% initialize vars and pars.
clear;
nx=500; nt=1800; %grid size and num time steps
v(1:nx)=0; u(1:nt,1:nx)=0; %initialise arrays
vm=v;sz=v;sx=v;mu=v;rho=v;ro1=v;muu=v;
tau_b=0.03; %source duration
nxs=60; %source position
nxr=200; %receiver position to record seismogram
dt=.002; dx=10.; %sampling step time & space
plo=0;
% stiffness and density values:
% * case (1), the string is homogeneous, use scalar values for parameters.
%
% * case (2), the sting is inhomogeneous; define parameters as a function of position and
% then compute the harmonic mean (this is done to match continuity conditionsacross
% interfaces).
%%%%%%
% case (1):
% mu=1e10; rho=2500.;
%%%%%%
% case (2):
mu(1:200)=1e10; rho(1:200)=2500.;
% contrast is at position 61*dx.
mu(201:nx)=2.5e10; rho(201:nx)=2777.;
% effective parameters (for inhomogeneous media)
muu=mu; ro1=rho;
muu(2:nx-1)=2./(1./mu(2:nx-1)+ 1./mu(3:nx)) ;
ro1(2:nx-1)=2./(1./rho(2:nx-1)+1./rho(3:nx)) ;
%%%%%%
% ***time loop start***
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %compute velocities update
 ??
 % update displacement
 ??
 % compute stress update
 ??
 % implement boundary conditions
 ??
 % apply source: (this is velocity imposed at point nxs, Gaussian in time)
 ??
 % save for final output seismogram
 ??
 % plot each 50 iterations
 ??
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ***time loop end***
%
% plot time/space full simulation
??
%end

Part3: 2D wave

Follow the same steps as for the 1D wave equation, except here use one update for velocity, and one update for each component of the stress and . Namely:

3.11.) Write the 2D wave equation as a coupled of equations for velocity and stress rate (see slides)

3.12) Build space-time grid and discrete form.

3.13) Use the following dimensions and parameters:

m
space and time number of points: nx=500;nz=100; nt=1800;
stiffness and density values for : ; kg/m3.
stiffness and density values for : ; kg/m3.

3.14) Apply B.C.s with free surface (stress=0) at the top . Which components of the stress are zero at the free top surface? Is it , , or both? What is the implicit B.C. at the side and bottom boundaries?

3.15) Plot the wave propagation and the ‘seismograms’ (in this case more precisely velocigrams) monitored along the surface as a function of position. What is the shape of the reflected phase? Identify first arrivals: is it always the same phase which arrives first, independent of position? Can you identify the critically refracted arrival?

NOTE: The refracted wave propagates faster and overtakes the direct arrival after a critical distance. Critical refraction takes place when part of the wave propagates parallel to the interface between two different layers. In terms of rays this will happen when , where is the angle between the ray and the vertical in the top layer, and are the velocities of the waves in top and bottom layers respectively. As it continues to propagate along the interface, part of the refracted wave leaks back into the upper layer creating a planar upward propagating front.

Template

This code simulates wave propagation (SH or acoustic) in 2D. An interface with impedance contrast is introduced to create both refracted and reflected waves.

% initialize variables and and parameters
clear;nx=500;nz=100; nt=1800; ntr=(nt/10); %grid size and num time steps
v(1:nz,1:nx)=0; u(1:nt,1:nx)=0; %initialise arrays
vm=v;sz=v;sx=v;mu=v;rho=v;ro1=v;ro2=v;muu=v;
tau_b=0.02; %source duration
nxs=160; nzs=2; %source position horizontal and depth
nxr=200; nzr=2; %receiver position depth 0.5
dt=.001; dx=10.; %sampling step time & space
cc=0.93; % attenuation factor for sponge boundaries
plo=0;
nzi=40; %depth of interface
%
% assign stiffness density values top layer and bottom layer
 mu(1:nzi,:)=1e10; rho(1:nzi,:)=2500.;
% interface is at depth 61*dx; assign bottom layer
 mu(nzi+1:nz,:)=2.5e10; rho(nzi+1:nz,:)=2777.;
%
%effective parameters for inhomogeneous media (harmonic mean)
muu=mu; ro1=rho; ro2=rho;%
muu(2:nz,1:nx-1)=4./(1./mu(2:nz,1:nx-1)+ 1./mu(2:nz,2:nx)+ 1./mu(1:nz-1,2:nx)+ 1./mu(1:nz-1,1:nx-1)) ;
ro1(2:nz,1:nx-1)=2./((1./rho(2:nz,1:nx-1)+1./rho(2:nz,2:nx))) ;
ro2(2:nz,1:nx-1)=2./((1./rho(2:nz,1:nx-1)+1./rho(1:nz-1,1:nx-1))) ;
%
% ***time loop start***
% compute velocities by time extrapolation :
v(1:nz-1,1:nx-1)= ??
% -- BEGIN SPONGE BOUNDARIES --
% ADD THIS SECTION IN ONLY ONCE THE REST WORKS!...:
atf=0.07;nbo=70;wi=nbo/2.5;
for j=1:nz;
 for i=1:nbo;
 v(j,i)=v(j,i)*(1.-atf*exp(-((i-1)/wi)^2));
 sx(j,i)=sx(j,i)*(1.-atf*exp(-((i-1)/wi)^2));
 sz(j,i)=sz(j,i)*(1.-atf*exp(-((i-1)/wi)^2));
 end;
 for i=nx-nbo+1:nx;v(j,i)=v(j,i)*(1.-atf*exp(-((i-nx)/wi)^2));
 sx(j,i)=sx(j,i)*(1.-atf*exp(-((i-nx)/wi)^2));
 sz(j,i)=sz(j,i)*(1.-atf*exp(-((i-nx)/wi)^2));
 end;
end
for j=nz-nbo+1:nz;
 for i=1:nx;
 v(j,i)=v(j,i)*(1.-atf*exp(-((j-nz)/wi)^2));
 sx(j,i)=sx(j,i)*(1.-atf*exp(-((j-nz)/wi)^2));
 sz(j,i)=sz(j,i)*(1.-atf*exp(-((j-nz)/wi)^2));
 end
end;
% -- END SPONGE BOUNDARIES--
%
% compute stress by interpolation
sx(1:nz,2:nx)= ??
sz(2:nz,1:nx)= ??
% apply free traction condition
??
% Source: (this is force in direction z, gaussian in time)
rf=0. ;
rf=(1./tau_b)*exp(-((it*dt-2.*tau_b)/tau_b)^2) ;
v(nzs,nxs)=v(nzs,nxs)+rf ;
% save output seismogram
u(it,:)=v(2,:);
if (plo < 100);
 plo=plo+1
else
 plo=0;
end
figure (1);
vp=v; %this is a dummy copy of the velocity, used for plotting only.
% Change it here to show :
vp(nzi:nzi+1,:)=100; % the interface position
vp(nzs:nzs+2,nxs:nxs+2)=100; % the source position
??
colormap(hsv);set(gca,'DataAspectRatio',[1 1 1]);
legend('v','Location','South');% colorbar;
set(gca,'DataAspectRatio',[1 1 1]); axis tight;
%
figure(2);
??
view([0 0 -1]); caxis manual;caxis ([0 2]);
colormap(flipud(gray));set(gca,'DataAspectRatio',[1 1 1]);
%legend('u','Location','South');% colorbar;
set(gca,'DataAspectRatio',[1 1 1]); axis square;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ***End time loop ***

Part 4 (advanced): higher order interpolation

3.16) Write series expansion of v at x up to second order

3.17) Write series expansion of v at x+dx up to second order

3.18) Compute the difference

3.19) Build space-time grid and discrete form

3.20) Program into matlab and plot

Further notes on higher order interpolation:

Taylor series, second and fourth order: Use the Taylor series expansion developed up to the second term to show that the ‘centred’ derivative of a function f(x) can be approximated as:

(In the limit h 0, this is the definition of the derivative). Use the same approximation again to show that:

and (f’ denotes first derivative, f” second and f”’ third derivative w/r to x):

Develop Taylor series up to the fourth term and use the above approximations to show that the ‘centred’ derivative of a function f(x) can be also approximated as:

Now assuming that the step h =dx/2, we obtain the definition of “second order” and “fourth order” discrete derivatives:

Application 1: the finite difference 1D wave equation:

The 1D wave equation may be written in terms of velocity and stress :

where s(.) is a forcing source term, and we arbitrarily chose the position and time . Using the second order centred approximation of and we obtain:

Using Hooke’s law time derivative , and using the fourth order approximation in space we obtain for :

Published with MATLAB® R2016a

NNESMO Practical 3: The wave equation

Contents

Introduction

Part 1: Harmonic oscillator

Part 2: 1D wave on a spring

Part3: 2D wave

Part 4 (advanced): higher order interpolation

Further notes on higher order interpolation: