<p>The onset of cavitation in hydraulic devices is caused by rapid pressure variations leading to the formation and collapse of vapor bubbles which can occasionally provoke severe damage. Shape optimizations for cavitation reduction, based on single-phase flow models, prevent static pressure from dropping below the vaporization pressure, but they do not account for the influence of vapor in the flow. In this article, a continuous adjoint-based optimization method for two-phase turbulent cavitating flows, using the Volume of Fluid method, is developed within an in-house GPU-enabled CFD solver. The adjoint accommodates new terms arising from the differentiation of the source terms modeling cavitation, the liquid–vapor mixture properties, as well as the turbulence model. The sensitivity derivatives are expressed in terms of surface integrals, using an adjoint formulation which, prior to this work, has been applied only to single-phase flows. The comparison of the sensitivity derivatives computed by the continuous adjoint method and finite differences shows that the differentiation of the turbulence model equations, frequently omitted in two-phase flows, is necessary to predict accurate gradients. After the validation of both the primal and adjoint solvers, gradient-based constrained shape optimizations are performed in three cavitation-dominated flows around a 2D isolated hydrofoil and a 3D hemispherical head body, as well as inside a 2D hydraulic poppet valve.</p>