The harnessing of mechanical power from supersonic flows is constrained by physical limitations and substantial aerodynamic losses. Bladeless axial turbines are a viable alternative to extract power in such harsh conditions without restricting the operating conditions. In this paper, we present a shape optimization of the wavy surface of bladeless turbines to maximize the power extraction, while minimizing convective heat fluxes and pressure losses. First, a baseline geometry was defined and an experimental campaign was carried out on the baseline wavy surface of the bladeless turbine at supersonic conditions in the Purdue Experimental Aerothermal Lab. Pressure, heat flux and skin friction measurements were compared with the Reynolds Averaged Navier Stokes results. Afterwards, the evaluation routine which consisted of the blade generation, grid generation, solving, and post-processing was implemented within an evolutionary optimizer with a multi-objective function to maximize the pressure force and minimize heat flux and pressure loss. Finally, a three-dimensional assessment in terms of power, heat load and pressure drop was performed for the best performing geometry with the commercial solver CFD++ of Metacomp. Turbulence closure was provided with the k-omega-SST turbulence model. The annular chamber of the bladeless turbine consisted of an unstructured mesh of approx. 8–10 million grid points.