The theme of the research is on the use of the hypersingular boundary element method for the modelling of surface waves. Surface waves in solids are known to be partially reflected & transmitted and mode converted into body waves at stress discontinuities, which suggests that a formulation continuous in stress and strain might prove beneficial for modelling purposes. Such continuity can be achieved with a subparametric approach where the geometry is approximated using linear elements and the field variables, displacement and traction, are approximated using cubic Hermitian and linear shape functions respectively. The higher order polynomial for approximating displacement is intended to be a more accurate representation of the physics relating to surface wave phenomena, especially at corners, and thus, is expected to capture this behaviour with greater accuracy than the standard isoparametric approach. The subparametric approach affords itself to continuity in stress and strain by imposing a smoothness in the elements, which is not available to the isoparametric approach. As the attention is focused primarily on the modelling of surface waves on the boundary of a medium rather than the interior, the boundary element method lends itself appropriately to this end.A 2D semi analytical integration scheme is employed to evaluate the integrals appearing in the hypersingular boundary integral formulation. The integration scheme is designed to reduce the errors incurred when integrals with singular integrands are evaluated numerically. The scheme involves the application of Taylor expansions to formulate the integrals into two parts. One part is regular and is evaluated numerically and the other part is singular but sufficiently simple to be evaluated analytically. The scheme makes use of the aforementioned subparametric approach and is applied to linear elements for the use in steady state elastodynamic boundary element method problems. The steady state problem is used as it is a simplified problem and is sufficient to permit the investigation of surface vibration at a constant motion. The 2D semi analytical integration scheme presented can be naturally extended to 3D.A particular focus and novelty of the work is the application of different limiting approaches to determine the free terms common to boundary integral methods. The accurate numerical solution of hypersingular boundary integral equations necessitates the precise evaluation of free terms, which are required to counter discontinuous and often unbounded behaviour of hypersingular integrals at a boundary. The common approach for the evaluation of free terms involves integration over a portion of a circular/spherical shaped surface centred at a singularity and allowing the radius of the circle/sphere to tend to zero. This approach is revisited in order to ascertain whether incorrect results are possible as a consequence of shape dependency, which is a recognised issue for hypersingular integrals.Two alternative methods, which are shape invariant, are proposed and investigated for the determination of free terms. The first approach, the point limiting method, involves moving a singularity towards a shrinking integration domain at a faster rate than the domain shrinks. Issues surrounding the choice of approach, shrinkage rates and path dependency are examined. A related and second approach, the boundary limiting method, involves moving an invariant, but shrinking, boundary toward the singularity, again at a faster rate than the shrinkage of the domain. The latter method can be viewed as a vanishing exclusion zone approach but the actual boundary shape is used for the boundary of the exclusion zone. Both these methods are shown to provide consistent answers and can be shown to be directly related to the result obtained by moving a singularity towards a boundary, that is, by comparison with the direct method. Unlike the circular/spherical approach the two methods involve integration over the actual boundary shape and consequently shape dependency is not an issue. A particular highlight of the point limiting approach is the ability to obtain free terms in mixed formulation, which is not available to the circular/spherical approach.There are three numerical problems considered in this research. The first problem considers the longitudinal vibration of a square plate. This is a problem for which a known analytical solution exists and is used to verify the equation formulation and integration scheme adopted for the isoparametric and subparametric formulations. Both formulations are as accurate as each other and produce results that are in keeping with the analytical solution, thus instilling confidence in their predictions.The second problem considers the simulation of surface waves on a square plate. Various boundaries of a square plate have displacement conditions imposed on them as a result of surface wave propagation. The results indicate that the surface wave behaviour is not captured. However, the analytical solution does not make any consideration for the effects from corners; the analytical solution is for a Rayleigh wave propagating upon a planar surface. It does not take into account the wave phenomena encountered at corners. Therefore, these results cannot be used to validate the predictions obtained on the boundary of the problem considered. The purpose of this problem is to illustrate the impact of corners on the surface wave propagation. Sensitivity studies are conducted to illustrate the effect of corners on the computed solution at the boundary.The final problem considers the simulation of surface waves on a circular plate. Various portions of the boundary of the circular plate have displacement conditions imposed on them as a result of surface wave propagation on curved surfaces. The results indicate that the isoparametric and subparametric predictions are similar to one another. However, both displacement profiles predict the presence of other waves. Given the multi faceted nature of the mesh, the computed solution is picking up mode conversion and partial reflection & transmission of surface waves. In reality, this is not expected as the surface of the boundary is smooth. However, due to the discretisation there are many corners in this problem.