We investigate properties of exponential operators preserving the particle number using combinatorial methods developed in order to solve the boson normal ordering problem. In particular, we apply generalized Dobinski relations and methods of multivariate Bell polynomials which enable us to understand the meaning of perturbation-like expansions of exponential operators. Such expansions, obtained as formal power series, are everywhere divergent but the Pade summation method is shown to give results which very well agree with exact solutions got for simplified quantum models of the one mode bosonic systems.