Moving contact line problem plays an important role in fluid-fluid interface motion on solid surfaces. The problem can be described by a phase-field model consisting of the coupled Cahn-Hilliard and Navier-Stokes equations with the generalized Navier boundary condition (GNBC). In this thesis, we generalize the GNBC to surfaces with complex geometry and introduce a finite element method on unstructured 3D meshes. Two efficient numerical methods are presented, including a linear decoupled scheme and a linearized coupled scheme. Numerical experiments are carried out to validate the effectiveness and efficiency of the proposed schemes. Accurate simulation of the interface and contact line motion requires very fine meshes, and the computation in 3D is even more challenging. Thus, the use of high performance computers and scalable parallel algorithms are indispensable. A highly parallel solution strategy using different solvers for different components of the discretization is presented. Parallel performances show that the strategy is scalable for 3D problems on a supercomputer with a large number of processors. We apply the proposed schemes and solution algorithms to study three important applications, particularly for those phenomena that can not be achieved by 2D simulations: 1) a droplet spreading on a rough surface; 2) flow features for a solid object impacting on a liquid surface; 3) the capillary force hysteresis and pinning-depinning events by direct simulations of a fiber intersecting a liquid-air interface.