Vertical surface piercing cylinders, such as typical coastal wind turbine foundations and basic elements of many coastal structures, are often exposed to combined loading from waves and currents. Accurate prediction of hydrodynamic loads on a vertical cylinder in a combined wave-current flow is a challenging task. This work describes and compares two different approaches for numerical modelling of the interaction between focussed wave groups and a sheared current, and then their interactions with a vertical piercing cylinder. Both approaches employ an empirical methodology to generate a wave focussed at the location of the structure in the presence of sheared currents and use OpenFOAM, an open source Computational Fluid Dynamics (CFD) package. In the first approach, the empirical wave-on-current focussing methodology is applied directly in the OpenFOAM domain, replicating the physical wave-current flume. This approach is referred to as the Direct Method. In the second approach, a novel Lagrangian model is used to calculate the free surface elevation and flow kinematics, which are then used as boundary conditions for a smaller 3-D OpenFOAM domain with shorter simulation time. This approach is referred to as the Coupling Method. The capabilities of the two numerical methods have been validated by comparing with the experimental measurements collected in a wave-current flume at UCL. The performance of both approaches is evaluated in terms of accuracy and computational effort required. It is shown that both approaches provide satisfactory predictions in terms of local free surface elevation and nonlinear wave loading on the vertical cylinders with an acceptable level of computational cost. The Coupling Method is more efficient because of the use of a smaller computational domain and the application of the iterative wave-current generation in the faster Lagrangian model. Additionally, it is shown that a Stokes-type perturbation expansion can be generalized to approximate cylinder loads arising from wave groups on following and adverse sheared currents, allowing estimation of the higher-order harmonic shapes and time histories from knowledge of the linear components alone.