In this paper, we present an adaptive moving mesh algorithm for meshes of unstructured polyhedra in three space dimensions. The algorithm automatically adjusts the size of the elements with time and position in the physical domain to resolve the relevant scales in multiscale physical systems while minimizing computational costs. The algorithm is a generalization of the moving mesh methods based on harmonic mappings developed by Li et al. [J. Comput. Phys., 170 (2001), pp. 562-588, and 177 (2002), pp. 365-393]. To make 3D moving mesh simulations possible, the key is to develop an efficient mesh redistribution procedure so that this part will cost as little as possible comparing with the solution evolution part. Since the mesh redistribution procedure normally requires to solve large size matrix equations, we will describe a procedure to decouple the matrix equation to a much simpler block-tridiagonal type which can be efficiently solved by a particularly designed multi-grid method. To demonstrate the performance of the proposed 3D moving mesh strategy, the algorithm is implemented in finite element simulations of fluid-fluid interface interactions in multiphase flows. To demonstrate the main ideas, we consider the formation of drops by using an energetic variational phase field model which describes the motion of mixtures of two incompressible fluids. Numerical results on two- and three-dimensional simulations will be presented.