Predictive simulation of two-phase flows with moving contact lines is challenging due to their inherent multi-physics and multi-scale nature. Directly simulating such flows incurs an enormous computational cost due to the widely disparate scales at the contact line. Moreover, simulations of viscosity-dominated two-phase flows with moving contact lines are often reported to be mesh-dependent due to the diverging viscous stress at the contact point. This dissertation addresses the above simulation issues and by building a numerical framework to enable large-scale 3D simulations of two-phase flows in complex geometries. By analyzing the weak form of the Navier-Stokes equations for a control volume adjacent to a wall with moving contact line, two unclosed terms are identified: a sub-grid scale (SGS) surface tension force and an SGS viscous force. A closure for the SGS surface tension force is first proposed and tested in a numerical framework for simulating two-phase flows with contact lines. This framework combines a conservative level set method to capture the interface and a conservative cut-cell immersed boundary method to handle complex geometries. Detailed verification tests confirm that simulations using this framework are discretely conservative, accurate, and robust. Secondly, a physics-based closure is derived for the SGS viscous force. Simulations these two SGS models are verified to be mesh-independent and physically accurate across a number of viscosity-dominated two-phase flows, including drop spreading on a horizontal plane and drop sliding down an inclined plane. Finally, the present approach is applied in the study of drop-fiber interactions and jet-wall interactions.