1.. -*- rst -*- -*- restructuredtext -*- 2 3.. This file should be written using the restructure text 4.. conventions. It will be displayed on the bitbucket source page and 5.. serves as the documentation of the directory. 6 7.. default-role:: math 8 9######################################################## 10PeWe - Particular solutions to the elastic wave equation 11######################################################## 12 13This software provides Matlab and Fortran routines that compute 14particular solutions to the elastic wave equation. 15 16Currently we provide solutions for the following problems: 17 18 1. Scattering by a circular cylindrical cavity in an elastic medium, 19 2. Scattering and diffraction by a circular cylindrical inclusion in an elastic medium, 20 3. Surface waves on the convex boundary of an elastic circular cylinder, 21 4. Surface waves on the concave boundary of a circular cylindrical cavity in an elastic medium. 22 23Below we describe how to use the provided routines and briefly explain 24the solutions, more comprehensive derivation of the solutions can be 25found in the paper `Formulae and Software for Particular Solutions to the Elastic Wave Equation in Curved Geometries`__. 26 27 28The Elastic Wave Equation in Cylindrical Geometries 29--------------------------------------------------- 30 31The radius of the cylinder is `a>0` and its axis is parallel to one of 32the coordinate axes, say :math:`z`. 33We consider the case when the material properties of the cylinder and 34the surrounding medium are different and we allow for either to be a 35vacuum. 36The displacement in the direction of the cylindrical axis is omitted, 37the remaining radial and azimuthal components of the displacement, `p(r,\theta,t)` and `q(r,\theta,t)`, are functions of the cylindrical coordinates, `r` and `\theta` and can be expressed as 38 39.. math:: 40 41 p = \phi_r + \frac{1}{r} \psi_\theta,\\ 42 q = \frac{1}{r}\phi_\theta - \psi_r, 43 44The functions `\phi` and `\psi` represent pressure and shear waves 45propagating with phase velocities `c_p` and `c_s`, which can be 46expressed in terms of the density, `\rho` of the elastic medium, 47Lame\'s first parameter, `\lambda` and the shear modulus `\mu` 48 49.. math:: 50 51 c_p = \sqrt{\frac{\lambda + 2 \mu}{\rho}}, c_s = \sqrt{\frac{\mu}{\rho}}, 52 53In Cartesian coordinates `x = r \cos(\theta)` and `y = r \sin(\theta)` and the horizontal and vertical displacement components `u(x,y,t)` and `v(x,y,t)` are given by 54 55.. math:: 56 57 u(x,y,t) = \cos(\theta) p(r,\theta,t) - \sin(\theta) q(r,\theta,t),\\ 58 v(x,y,t) = \sin(\theta) p(r,\theta,t) + \cos(\theta) q(r,\theta,t). 59 60These displacements are output by the routines outlined below. For 61examples of how the routines are called within Matlab see 62``usage_example.m`` and for Fortran, see the files starting with 63``main_`` in the fortran directory. 64 65Scattering by a circular cylindrical cavity in an elastic medium 66---------------------------------------------------------------- 67Particular solutions that represent a plane wave impinging on a cylidrical cavity and the resulting scattered wave fields are given by 68 69.. math:: 70 71 p(r,\theta,t) = e^{i\omega t}\sum_{n=0}^{\infty} \left(\phi_0 \epsilon_n (-i)^n \frac{d}{dr} J_n(\gamma_p r)+ A_n \frac{d}{dr} H_n^{(2)} (\gamma_p r) + B_n\frac{n}{r} H_n^{(2)} (\gamma_s r) \right) \cos(n \theta),\\ 72 q(r,\theta,t) = e^{i\omega t}\sum_{n=0}^\infty \left(\frac{-n}{r}\phi_0 \epsilon_n (-i)^n J_n(\gamma_p r) + A_n\frac{-n}{r} H_n^{(2)}(\gamma_p r)-B_n \frac{d}{dr} H_n^{(2)} (\gamma_s r)\right) \sin(n\theta). 73 74Here `J_n` and `H^{(2)}_n` are the Bessel function of order `n` and the Hankel function of the second kind of order `n`, respectively. The amplitude of the incoming plane wave is `\phi_0` and its temporal frequency is `\omega`. The remaining quantities `\gamma_p,\gamma_s, \epsilon_n,A_n,B_n` are determined so that the solution satisfies the elastic wave equation and the boundary conditions at the cylindrical cavity, for details se the references below. The routine ``[du, dv] = cylindrical_cavity(X,Y,t,omega,lambda,mu,rho)`` returns the real part of the resulting vertical and horizontal displacements in Cartesian coordinates. The routine is called in the following way 75 76.. code-block:: 77 78 function [du, dv] = cylindrical_cavity(X,Y,t,omega,lambda,mu,rho) 79 subroutine cylindrical_cavity(du,dv,X,Y,t,omega,lambda,mu,rho) 80 81Here ``X,Y`` are the location(s) where the solution is to be computed. 82In the Matlab implementation ``X,Y`` can be either scalars or ``nx`` 83by ``ny`` matrices, in the Fortran implementation they are assumed to 84be scalars. 85 86Scattering and diffraction by a circular cylindrical inclusion in an elastic medium 87----------------------------------------------------------------------------------- 88Particular solutions that represent the resulting diffacted wave field inside a cylindrical inclusion resulting from the incoming plane wave described above are given by 89 90.. math:: 91 92 p'(r,\theta,t) = e^{i\omega t}\sum_{n=0}^{\infty} \left(C_n \frac{d}{dr} J_n (\gamma_p' r) + D_n\frac{n}{r} J_n (\gamma_s' r) \right) \cos(n \theta),\\ 93 q'(r,\theta,t) = e^{i\omega t}\sum_{n=0}^\infty \left(C_n\frac{-n}{r} H_n^{(2)}(\gamma_p' r)-D_n \frac{d}{dr} J_n(\gamma_s' r)\right) \sin(n\theta). 94 95The quantities `C_n,D_n,\gamma_P',\gamma_S'` are determined by the material parameters and so that both `p,q` above and `p',q'` satisfies the conditions at the interface between the cylinrdical inclusion and the surrounding medium, for details see the references below. The routine ``[du, dv,du_p, dv_p] = cylindrical_inclusion(X,Y,X_p,Y_p,t,omega,lambda,mu,rho,lambda_p,mu_p,rho_p)`` returns the real part of the resulting vertical and horizontal displacements both outside and inside the cylindrical inclusion and in Cartesian coordinates. The routine is called in the following way 96 97.. code-block:: 98 99 function [du, dv,du_p,dv_p] = cylindrical_inclusion(X,Y,X_p,Y_p,t,omega,lambda,mu,rho) 100 subroutine cylindrical_inclusion(du,dv,du_p,dv_p,X,Y,t,omega,lambda,mu,rho) 101 102Here ``X,Y`` and ``X_p,Y_p`` are the location(s) outside and inside 103the cylindrical inclusion, respectively, where the solution is to be 104computed. In the Matlab implementation both ``X,Y`` and ``X_p,Y_p`` 105can be either scalars or ``nx`` by ``ny`` matrices. The Fortran 106implementation requires the coordinates to be scalars. 107 108Surface waves on the convex boundary of an elastic circular cylinder 109-------------------------------------------------------------------- 110 111Particular solutions that represent time - harmonic circumferential 112waves that travel in the tangential direction along boundary of the 113cylinder with phase velocity `c` are (expressed as radial and azimuthal displacements) 114 115.. math:: 116 117 p(r,\theta,t) = \left(A \frac{d J_{ka}(c \eta_p r)}{dr} + B \frac{i ka}{r} J_{ka}(c\eta_s r)\right) e^{i ka(c t + \theta)},\\ 118 q(r,\theta,t) = \left(A\frac{ik a}{r} J_{ka}(c \eta_p r) - B \frac{dJ_{ka}(c\eta_s r)}{dr} \right) e^{i ka(c t + \theta)}. 119 120The routine ``surface_waves_convex(X,Y,t,n,c,B)`` returns these 121waves. The routine is called in the following way 122 123 124.. code-block:: 125 126 function [du, dv] = surface_waves_convex(X,Y,t,n,c,B) 127 subroutine surface_waves_convex(du,dv,X,Y,nx,ny,t,n,c,B) 128 129 130Here ``X,Y`` are the location(s) where the solution is to be 131computed. In the Matlab implementation ``X,Y`` can be either scalars 132or ``nx`` by ``ny`` matrices. In the Fortran implementation they are 133required to be scalars. 134 135 136Surface waves on the concave boundary of a circular cylindrical cavity in an elastic medium 137------------------------------------------------------------------------------------------- 138 139 140Particular solutions that represent time - harmonic circumferential 141waves that travel in the tangential direction along boundary of the 142cylinder with phase velocity `c` are (expressed as radial and azimuthal displacements) 143 144.. math:: 145 146 p(r,\theta,t) = \left(A \frac{d H^{(2)}_{ka}(c \eta_p r)}{dr} + B \frac{i ka}{r} H^{(2)}_{ka}(c\eta_s r)\right) e^{i ka(c t + \theta)},\\ 147 q(r,\theta,t) = \left(A\frac{ik a}{r} H^{(2)}_{ka}(c \eta_p r) - B \frac{dH^{(2)}_{ka}(c\eta_s r)}{dr} \right) e^{i ka(c t + \theta)}. 148 149The routine ``surface_waves_concave(X,Y,t,n,c,B)`` returns these 150waves. The routine is called in the following way 151 152 153.. code-block:: 154 155 function [du, dv] = surface_waves_concave(X,Y,t,n,c,B) 156 subroutine surface_waves_concave(du,dv,X,Y,nx,ny,t,n,c,B) 157 158 159Here ``X,Y`` are the location(s) where the solution is to be 160computed. In the Matlab implementation ``X,Y`` can be either scalars 161or ``nx`` by ``ny`` matrices. In the Fortran implementation they are 162required to be scalars. 163 164For this solution the Fortran routines makes use of the subroutine 165``ZBESH`` which can be obtained from Netlib or from 166``https://github.com/JuliaLang/openspecfun.git``. 167 168 169Licence 170------- 171PeWe is written by Kristoffer Virta, and Daniel Appelo 172and released under the GNU General Public Licence version 3 (or later). 173 174__ https://bitbucket.org/appelo/pewe/raw/51b7c354bd5120686a66a4bae8b63086feb01be4/documentation/part_sol.pdf 175