Conventional finite element methods (FEM) have been used for many years for the solution of harmonic wave problems. To ensure accurate simulation, each wavelength is discretised into around ten nodal points, with the finite element mesh being updated for each frequency to maintain adequate resolution of the wave pattern. This technique works well when the wavelength is long or the model domain is small. However, when the converse applies and the wavelength is small or the domain of interest is large, the finite element mesh requires a large number of elements, and the procedure becomes computationally expensive and impractical. The principal objective of this work is to accurately model two-dimensional Helmholtz problems with the Partition of Unity Finite Element Method (PUFEM). This will be achieved by applying the plane wave basis decomposition to the wave field. These elements allow us to relax the traditional requirement of around ten nodal points per wavelength and therefore solve Helmholtz wave problems without refining the mesh of the computational domain at each frequency. Various numerical aspects affecting the efficiency of the PUFEM are analysed in order to improve its potential. The accuracy and effectiveness of the method are investigated by comparing solutions for selected problems with available analytical solutions or to high resolution numerical solutions using conventional finite elements. First, the use of plane waves or cylindrical waves in the enrichment process is assessed for wave scattering problems involving a rigid circular cylinder in both near field and far field. In the far field, the cylindrical waves proved to be more effective in reducing the computational effort. But given that the plane waves are simpler to analytically integrate for straight edge elements, during the finite element assembling process, they are retained for the remaining of the thesis. The analysed numerical aspects, which may affect the PUFEM performance, include the conjugated and unconjugated weighting, the geometry description, the use of non-reflecting boundary conditions, and the h-, p- and q-convergence. To speed up the element assembling process at high wave numbers, an exact integration procedure is implemented. The PUFEM is also assessed on multiple scattering problems, involving sets of circular cylinders, and on exterior wave problems presenting singularities in the geometry of the scatterer. Large and small elements, in comparison to the wavelength, are used with both constant and variable numbers of enriching plane waves. Last, the PUFEM resulting system is iteratively solved by using an incomplete lower and upper based preconditioner. To further enhance the efficiency of the iterative solution, the resulting system is solved into the wavelet domain. Overall, compared to the FEM, the PUFEM leads to drastic reduction of the total number of degrees of freedom required to solve a wave problem. It also leads to very good performance when large elements, compared to the wavelength, are used with high numbers of enriching plane waves, rather than small elements with low numbers of plane waves. Due to geometry detail description, it is practical to use both large and small elements. In this case, to keep the conditioning within acceptable limits it is necessary to vary the number of enriching plane waves with the element size.