Abstract The use of mesh refinement techniques is becoming more and more popular in computational fluid dynamics, from multilevel approaches to adaptive mesh refinement. In this paper we present a new method to interpolate the coarse velocity field which is based on an optimal approach and is characterized by a constrained minimization of an objective functional. The functional contains the sum of the square difference between the velocity components and their target average value subject to a number of divergence-free constraints. In this work we describe this approach in two- and three-dimensional geometries with different discrete velocity field configurations. This technique is applied to a multilevel Volume-of-Fluid (VOF) method where the volume fraction function is used to reconstruct and advect the interface between two immiscible phases. The coarse velocity field is interpolated to a fixed fine grid with the optimal approach over a given number of refinement levels. The results of several kinematic tests are presented, where the mass and geometrical errors are compared with those obtained with refined velocity fields interpolated with a simple midpoint rule.