Using a sampling-based reliability-based design optimization method, we present a shape reliability-based design optimization method for coupled fluid–solid interaction problems. For the fluid–solid interaction problem in arbitrary Lagrangian–Eulerian formulation, a coupled variational equation is derived from a steady state Navier–Stokes equation for incompressible flows, an equilibrium equation for geometrically nonlinear solids, and a traction continuity condition at interfaces. The fluid–solid interaction problem is solved using the finite element method and the Newton–Raphson scheme. For the fluid mesh movement, we formulated and solved a pseudo-structural sub-problem. The shape of the solid is modeled using the Non-Uniform Rational B-Spline (NURBS) surface, and the coordinate components of the control points are selected as random design variables. The sensitivity of the probabilistic constraint is calculated using the first-order score functions obtained from the input distributions and from the Monte Carlo simulation on the surrogate model constructed by using the Dynamic Kriging method. The sequential quadratic programming algorithm is used for the optimization. In two numerical examples, the proposed optimization method is applied to the shape design problems of solid structure which is loaded by prescribed fluid flow, and this proves that the sampling-based reliability-based design optimization can be successfully utilized for obtaining a reliable optimum design in highly nonlinear multi-physics problems.