MATLAB Examples

The "dispersion relation" toolbox for surface waves

F. Moisy, 28 april 2010.

Laboratory FAST, University Paris Sud.


About the dispersion relation of surface waves

This toolbox simply provides an easy way to work with the dispersion relation of surface waves, given by

   omega(k) = sqrt ( tanh(k*h0) * (g*k + gamma*k^3/rho))

where omega is the pulsation (in rad/s), k the wavenumber (in 1/m), h0 the depth (in m), g the gravity (in m^2/s), gamma the surface tension (in N/m) and rho the fluid density (in kg/m^3).

The following example simply plots the wave frequency as a function of wavelength:

lambda = 0:0.1:10;  % wavelength (in m)
freq = omega(2*pi./lambda)/(2*pi);
xlabel('\lambda (m)');
ylabel('f (Hz)');

Inverting the dispersion relation

The function kfromw allows one to invert the dispersion relation, i.e. to give the value of omega for a given value of k.

For infinite depth, kfromw inverts the cubic polynomial. For finite depth, a zero-finding method is used, starting from the infinite depth solution.

freq = 10:50;  % in Hz;
lambda = 2*pi./kfromw(2*pi*freq);
plot(freq, 1000*lambda);
xlabel('f (Hz)');
ylabel('\lambda (mm)');

Setting the wave parameters

By default, the physical parameters (liquid densities, surface tension, etc.) are set for an air-water interface under usual conditions, with a water layer of infinite depth ("deep water waves"). The default settings are given by the function wave_parameter:

p = wave_parameter
p = 

        g: 9.8100
    gamma: 0.0720
      rho: 1000
       h0: Inf

The structure p can be passed as an input argument to omega and kfromw.

The default settings can be changed as follows: here the surface tension is set to 0.06 N/m and the depth to 0.1 m:

p = wave_parameter('gamma',0.05,'h0',0.1)
p = 

        g: 9.8100
    gamma: 0.0500
      rho: 1000
       h0: 0.1000

Re-do the same plot as above with the new settings:

lambda = 0:0.1:10;
freq = omega(2*pi./lambda, p)/(2*pi);     % note the second parameter p
xlabel('\lambda (m)');
ylabel('f (Hz)');

Phase and group velocity

The function omega also returns, as additional output arguments, the phase and group velocity of the wave:

k = linspace(1,1e4,1e4);
[w,c,cg] = omega(k);
lambda = 2*pi./k;
xlabel('\lambda, m')
ylabel('c, c_g,  m/s');
legend('Phase velocity c','Group velociy c_g');

Kelvin wake

A last function of the toolbox, kckg, returns the capillary wavenumber kc and gravity wavenumber kg for a given phase velocity U. This is useful in the problem of the Kelvin wake, in which the range of excited wavenumbers is given by equating the disturbance velocity and the phase velocity. Here again, the wave parameters stored in the structure p may be passed as an additional input argument to kckg.

v = linspace(0,1,1000);  % disturbance velocity, m/s
[kc, kg] = kckg(v);      % wavenumbers selected by the velocity
xlabel('disturbance velocity V, m/s');
ylabel('k_c,  k_g,  m^{-1}')
xlim([.1 1])
legend('Capillary branch','Gravity branch');

Rayleigh-Taylor instability

If the gravity is negative, omega may return a complex value. The imaginary part corresponds to the growthrate of the Rayleigh-Taylor instability (a denser fluid above a lighter fluid).

p = wave_parameter('g',-10);
k = linspace(100, 1000, 100);
s = omega(k,p);
lambda = 2*pi./k;
xlabel('\lambda (mm)');
legend('Growth rate (unstable)','Pulsation \omega (stable)');

Here the most unstable wavelength is k = 29 mm.


To learn more about this toolbox, see the help:

help dispersion_relation
  Dispersion Relation Toolbox (F. Moisy)
  Version 1.01 28-Apr-2010
  Dispersion relation
    omega          - Dispersion relation w(k) for surface waves
    kfromw         - Inverse dispersion relation k(w) for surface waves
    wave_parameter - Set the parameters for the dispersion relation
    kckg           - Wavenumbers selected by a given velocity
  Various plot examples
    plot_lambda    - Wavelength versus frequency
    plot_c         - Phase and group velocities versus wavenumber.
    plot_kckg      - Gravity and capillary wavenumbers of the Kelvin wake