In this work we investigate a phase field model for damage processes in two-dimensional viscoelastic media with non-homogeneous Neumann data describing external boundary forces. In the first part we establish global-in-time existence, uniqueness, a priori estimates and continuous dependence of strong solutions on the data. The main difficulty is caused by the irreversibility of the phase field variable which results in a constrained PDE system. In the last part we consider an optimal control problem where a cost functional penalizes maximal deviations from prescribed damage profiles. The goal is to minimize the cost functional with respect to exterior forces acting on the boundary which play the role of the control variable in the considered model. To this end, we prove existence of minimizers and study a family of "local" approximations via adapted cost functionals.