The main aim of this paper is to present new sets of non-perturbative fragmentation functions (FFs) for D0 and D+ mesons at next-to-leading (NLO) and, for the first time, at next-to-next-to-leading order (NNLO) in the MS¯¯¯¯¯¯¯ factorization scheme with five massless quark flavors. This new determination of FFs is based on the QCD fit to the OPAL experimental data for hadron production in the electron-positron single-inclusive annihilation (SIA). We discuss in detail the novel aspects of the methodology used in our analysis and the validity of obtained FFs by comparing with previous works in literature which have been carried out up to NLO accuracy. We will also incorporate the effect of charmed meson mass corrections into our QCD analysis and discuss the improvements upon inclusion of these effects. The uncertainties in the extracted FFs as well as in the corresponding observables are estimated using the “Hessian” approach. For a typical application, we use our new FFs to make theoretical predictions for the energy distributions of charmed mesons inclusively produced through the decay of unpolarized top quarks, to be measured at the CERN LHC. As a result of this analysis, suggestions are discussed for possible future studies on the current topic to consider any theory improvements and other available experimental observables.