Liquid metal (LM) plasma-facing components (PFCs) have the potential to alleviate some of the difficult constraints of solid PFCs. The solid substrate that supports these liquids would see only neutron irradiation and is protected from the plasma. Lithium or lithium eutectics would have a high affinity for tritium and deuterium at low operating temperatures and would provide a low recycling environment for the core plasma that seems to provide significant confinement improvements in existing plasma experiments. Liquid metals have sufficient thermal conductivity for good thermal conduction and can provide vapor shielding under transient heat loads. On the other hand, electrically conducting liquids will have magnetohydrodynamic (MHD) interactions that can disturb the surface allowing material to enter the plasma and can laminarize the flow, thereby reducing turbulence and convective heat transport. Sputtering and evaporation must be controlled by limiting the operating temperature. The substrate material must remain covered to survive, and the design for how the LM enters and exits the plasma chamber or divertor region needs to be identified. The present contribution introduces a liquid lithium free-surface cooling system into a future fusion device design developed for the Fusion Nuclear Science Facility (FNSF) program. The design includes LM PFCs covering convex divertor surfaces as well as lithium supply and removal lines. Three-dimensional analysis of the proposed design was performed using a highly customized version of the ANSYS CFX code applying methods of computational fluid dynamics. MHD was introduced using a magneto vector potential approach, allowing a natural interface. Free-surface flow was included using a volume of fluid approach with surface tension on the fluid vacuum interface. Electromagnetic equations were solved for LM as well as for solid components and vacuum. Special stabilization procedures were derived and applied to improve the convergence of the momentum equations with the source terms due to Lorentz force and surface tension. Conjugate heat transfer analysis was performed in LM and solid components. The numerical model was validated using analytical MHD solutions for high Hartmann number flow as well as relevant experimental results. Numerical solutions were obtained for the design configuration including FNSF geometry and Fusion Energy System Studies device parameters obtained from the system studies. Results of the analysis allow determination of the parameter range for the operation of the proposed design as well as design optimization.