Heat sources for ORC turbines typically feature variable energy sources such as WHR (Waste Heat Recovery) and solar energy. Advanced uncertainty quantification and robust optimization methodologies could be used during the ORC turbines design process in order to account for multiple uncertainties. This study presents an original robust shape optimization approach for ORC blade turbines, to overcome the limitation of a deterministic optimization that neglects the effect of uncertainties of operating conditions or design variables. Starting from a baseline blade, we search for an optimal shape that maximizes the 5% quantile of the expander isentropic efficiency, which is evaluated by means of an Euler 2D simulation. Real-gas effects are modeled through the use of a Peng-Robinson-Stryjek-Vera equation of state. The 5% quantile of the expander isentropic efficiency is estimated using a tail probability strategy: points are iteratively added on the failure branches in order to build a reliable metamodel from which a Monte-Carlo sampling method is used. In order to speed-up the optimization process, an additional Gaussian Process model is built to approximate the isentropic efficiency. The robustly optimized ORC turbine shape is finally compared to the initial configuration and the deterministic optimal shape.