In this study, the effect of the proposed binarization method on the topology optimization analysis of two-fluid flows with free surfaces was investigated. The oscillation of the liquid content of a vessel, referred to as sloshing, was numerically analyzed, and density-based topology optimization analysis was performed to suppress sloshing. In density-based topology optimization, the emergence of grayscale, which is neither a structural nor a fluid domain, is problematic. Therefore, in this study, the hyperbolic tangent function used in the sigmoid transformation to represent the structural and fluid domains was improved. In the early stages of the topology update, grayscale is allowed and the curves of the sigmoid transformation are loosely connected; however, as the number of iterations increases, the curves become steeper to eliminate grayscale. The optimization analysis was performed using a grayscale elimination method previously described in the research by Aage et al. (2008) and the proposed stepwise binarization method. The governing equations for the flow field consisted of the Navier-Stokes equations for unsteady incompressible viscous flow and the continuity equation. The interface was determined using the density function method. When dealing with incompressible viscous flow, caution should be exercised during the discretization process owing to the numerical instabilities caused by advection dominance and incompressibility conditions. To eliminate numerical instabilities, the governing equations were discretized based on the Lagrange-Galerkin and penalty methods. The numerical results confirmed that the grayscale was eliminated when the stepwise binarization method was employed.