Optimization with nonnegative orthogonality constraints has wide applications in machine learning and data sciences. It is NP-hard due to some combinatorial properties of the constraints. We first propose an equivalent optimization formulation with nonnegative and multiple spherical constraints and an additional single nonlinear constraint. Various constraint qualifications, the first- and second-order optimality conditions of the equivalent formulation are discussed. We design a class of exact penalty models in which the nonnegative and multiple spherical constraints are kept. The penalty models are exact if the penalty parameter is sufficient large other than going to infinity. A practical penalty algorithm with rounding technique is then developed. It uses a second-order method to approximately solve a series of subproblems with nonnegative and multiple spherical constraints. Extensive numerical results on the projection problem, orthogonal nonnegative matrix factorization problems and the K-indicators model show the effectiveness of our proposed approach.