It is well known that second order effects may in many cases be important for the nonlinear hydrodynamic problems arising in ocean engineering. Despite considerable efforts having been made in the past in calculating second order unsteady forces, similar studies are rare for the actual second order velocity potential itself, which is important for the understanding of wave kinematics. A mathematical model has been developed for the calculation of the second order sum frequency diffraction potential for fixed bodies in waves. It is believed that a first step towards the solution of the second order problem is the accurate evaluation of the first order quantities. By the use of Green's second identity, the first order problem can be cast into the form of a Fredholm integral equation and then solved by the Boundary Element Method. Some new developments based on this technique have been undertaken in this work, and as a result, there is a major improvement in the accuracy of the first order analysis. For the second order problem, the solution procedures are similar to those used for the first order problem except that special techniques have been developed to calculate efficiently the additional free surface integral which decays slowly to infinity in a highly oscillatory manner. In addition, an effective method has also been implemented to calculate the second derivative term in the free surface integral. From the numerical results presented, a number of interesting findings are illustrated. A closed form expression for a vertical circular cylinder has also been developed which not only furnishes a valuable check on the general numerical model but also provides some physical explanation for the second order phenomena. Moreover, it has been used to investigate some theoretical problems which (in the past) have caused confusion and error in the second order analysis. They are mainly associated with the troublesome nonhomogeneity presented in the free surface boundary condition.