Robust curve design on surface meshes with flexible controls is useful in a wide range of applications but remains challenging. Most existing methods fall into one of the two strategies: one is to discretize a curve into a polyline which is then optimized, and the other is to directly design smooth splines on meshes. While the former approach usually needs a sufficiently dense sampling of curve points which is computational costly, the latter approach relaxes the sampling requirement but suffers from the lack of user control. To tackle these problems, we proposed a variational method for designing feature-aware B-spline curves on surface meshes. Given the recent advances in shell space construction methods, we could relax the B-spline curve inside a simplified shell mesh and evaluate its distance to the surface using equipped bijective mapping. To effectively minimize the distance between the curve and the surface, with additional controls in the form of both internal and external constraints, we applied the interior point method and adaptively insert knots of the spline to increase its freedom and adjust the weighting during the iterations. When the curve is close enough to the surface, it can be efficiently sampled at any resolution and robustly projected to the surface. Experiments show that our method is more robust, has higher flexibility, and generates smoother results than existing methods.