An adjoint-based topology optimization method is developed for multiphase flows with phase change due to boiling for the applications such as optimization of the cooling efficiency of two-phase immersion cooling of batteries, power electronics and similar applications. Simulation of the two-phase flow is performed using a mixture model. Both the liquid and the vapor phases are assumed to be incompressible, whereas the mixture behaves as a compressible fluid; Laminar flow is assumed based on the Reynolds number present. The phase-change is modeled with Chen’s correlation, since nucleate boiling is expected as dominating boiling regime. For the optimization problem, a new formulation that considers the two-phase flow terms in the adjoint terms has been derived and implemented within the framework of the commercial CFD software AVL Fire. As objective functions, the overall pressure drop and the overall enthalpy difference through a flow channel where boiling occurs, are considered. For testing the adjoint method, a case study of a 2D heat sink is considered. The adjoint-based method for multiphase flows with phase change due to boiling is non-existent in the literature and is a novel contribution from the present work.